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ABSTRACT 


This  report  contains  a  study  of  joints  made  from  fiber  composite 
materials.  Several  types  of  joints  are  considered,  such  as  bolted 
and  adhesive  joints. 

Elements  for  a  rational  joint  design  are  presented,  as  well  as  test 
results.  A  finite  element  computer  program,  used  to  obtain  stress 
distributions  in  composite  joints,  is  included. 
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TECHNICAL  DISCUSSION 


BOLTED  JOINTS:  GENERAL  CONSIDERATIONS 

In  this  section,  the  fundamental  problems  of  the  design  of  bolted  joints 
in  composite  materials  (particularly,  glass-epoxy)  are  considered, 

A  bolted  joint,  Figure  1,  must  withstand,  in  general,  an  axial  force  P 
and  a  transverse  force  Q  .  Other  forces,  out  of  the  plane,  can  exist, 
but  in  most  cases  their  magnitudes  are  negligible  in  comparison  with  P 
and  Q  . 


Figure  l.  General  Bolted  Joint, 

Because  of  its  construction,  the  fiber  orientation  in  a  bolted  joint  is 
practically  confined  to  the  one  shown  in  Figure  1.  Fibers  through 
cross-section  A  are  continuous,  while  those  through  cross-section  B 
are  discontinuous  due  to  the  presence  of  the  bolt  hole.  When  the  axial 
load  P  is  tensile,  as  shown  in  Figure  1,  the  fibers  A  act  in  tension, 
but  the  fibers  B,  near  the  hole,  are  nearly  free  of  stresses.  Because  of 
this  stress  difference,  a  shear  stress  is  developed  along  the  fiber 
direction  starting  at  the  point  M  .  When  the  load  P  is  increased,  the 
shear  stress  at  M  produces  the  failure  of  the  matrix  of  the  composite. 

A  crack  is  propagated,  as  shown  in  Figure  2,  simultaneously  with  increas¬ 
ing  P  .  Test  results  indicate  that  the  magnitude  of  P  sufficient  to 
initiate  the  shear  crack  at  M  is  only  one-fourth  to  one-fifth  of  the 
value  that  produces  the  final  failure.  Once  the  shear  crack  is  completely 
developed,  the  joint  behaves  as  a  simple  loop  joint  for  tensile  load  P  , 
with  the  darkened  fibers  of  Figure  3  in  tension  and  the  central  region 
practically  free  of  stresses. 

If  the  joint  of  Figure  3  is  unloaded  and  then  reloaded  up  to  the  previous 
load  level,  no  further  changes  in  the  stress  distribution  are  observed. 
Thus,  the  shear  crack  has  no  further  effect  on  the  behavior  of  the  joint 
once  it  is  developed.  In  fact,  the  location  of  the  shear  crack  can  be 
controlled  if  the  fibers  A  are  not  bonded  to  the  fibers  B  during  the 
fabrication  process. 
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Figure  2.  Shear  Crack  in  Bolted  Joint 

After  First  Loading  in  Tension. 


Figure  3.  Joint  in  Tension  After  the  Propagation 
of  the  Shear  Crack. 
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If  the  tension  load  P  is  increased  after  the  propagation  of  the  shear 
crack,  a  final  failure  of  the  joint  will  appear  in  the  neighborhood  of 
the  point  N  (Figure  3).  The  stresses  at  the  point  N  are  as  shown  in 

Figure  4:  (a  tangential  tensile  stress  a  acting  along  the  fiber  direction 

0 

a 


Fiber  Orientation 


Figure  4.  Stresses  at  the  Points  N  (See  Figure  3). 


and  a  radial  compressive  stress  ,  perpendicular  to  the  fiber  direction). 

Because  the  stresses  and  a ^  act  simultaneously,  the  failure  of  the 

joint  is  actually  produced  by  combination  of  the  effects  of  each  stress* 

If  the  radius  a  of  the  hole  is  large  in  comparison  with  t  (see  Figure  5), 

the  stress  aa  is  much  larger  than  o  .  In  this  case,  the  failure  mode  of 
u  r 


Figure  5.  Nomenclature  for  the  Bolted  Joint, 
the  joint  is  a  tensile  failure  of  the  fibers,  as  represented  in  Figure  6. 
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If,  on  the  other  hand,  the  radius  a  is  much  smaller  than  t  ,  the 
radial  stress  ar  is  greater  than  qq  and  the  typical  failure  of  the 
joint  is  as  indicated  in  Figure  7.  In  this  case, 


Figure  7.  Typical  Failure  of  a  Joint  With  a«  t. 


the  failure  is  produced  essentially  by  the  compressive  radial  stress  . 
Only  the  region  near  the  bolt  is  affected  by  this  type  of  failure. 

From  the  previous  description  it  should  be  evident  that,  with  fixed 
external  radius  b  and  joint  thickness  h  ,  the  tensile  failure  load  P 
and  the  failure  mode  will  depend  upon  the  t/a  ratio.  Then,  P  will 
achieve  a  maximum  for  some  optimum  t/a  ratio.  By  performing  several 
tests  on  glass-epoxy  bolted  joints,  it  was  determined  that  this  optimum 
ratio  is  unity.  A  description  of  the  specimens  and  the  test  results  are 
given  in  the  next  section  of  this  report. 

Now,  let  us  consider  when  the  axial  load  P  is  compressive.  Even  for 
relatively  small  values  of  P,  the  failure  mode  depicted  in  Figure  8  was 
observed.  This  failure  is  due  to  a  wedging  effect  produced  by  the  bolt, 

Failure 


Figure  8.  Failure  Under  Compressive  Loading. 
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breaking  the  composite  in  the  direction  perpendicular  to  the  fibers.  In 
order  to  avoid  this  type  of  premature  failure  in  compression,  a  specially 

shaped  bushingfas  represented  in  Figure  9,  can  be  used. 


Bush ing 


Figure  9.  Joint  With  Special  Bushing  for  Compressive  Loading. 


When  the  compressive  load  P  is  applied  and  its  magnitude  increased,  a 
shear  crack  is  initiated  in  the  points  M  (see  Figure  9)  because  of  the 
stress  difference  between  the  compressed  fibers  of  the  central  region, 
shown  in  the  figure,  and  the  practically  stress-free  fibers  surrounding 
the  bushing.  After  the  completion  of  this  shear  crack  (Figure  9).  the 
fibers  of  the  central  region  behave  as  simple  compression  members.  The 
shape  of  the  contact  surface  between  the  bushing  and  the  compressed 
composite  avoids  any  wedging  effect  that  would  tend  to  separate  the  fibers. 
This  bushing  design  has  proven,  in  compression  tests,  to  be  very  effect¬ 
ive.  In  fact,  the  compression  ultimate  loading  in  joints  with  this  type 
of  bushing  is  always  greater  than  the  ultimate  tension  load. 

It  must  be  pointed  out  that  the  radius  a  is  the  external  radius  of  the 
bushing,  when  the  tension  design  is  employed. 

The  shear  force  Q  is,  in  most  cases,  only  a  fraction  of  P  .  The  joint 
design  previously  described  has  an  ultimate  shear  force  Q  of  about  one- 
tenth  of  the  ultimate  P  in  tension.  In  special  cases  when  Q  is,  for 
example, one  half  of  P  ,  either  the  joint  design  must  be  changed  or  the 
load  transmission  must  be  performed  in  another  manner  to  avoid  such  high 
shear  forces. 

The  ultimate  tensile  load  of  a  joint  with  a  ratio  t/a  =  1  has  been  found 
to  correspond  with  the  following  expression: 


P 


ult 


ult 


t  h 


a) 


where  h  is  the  joint  thickness. 
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This  formula  is  in  good  agreement  with  the  test  results.  It  should  be 
emphasized  that  it  is  generally  impossible  to  predict  the  ultimate  load 
of  a  composite  joint  from  a  stress  analysis  based  upon  a  linear  behavior 
of  the  material.  Such  analysis  must  be  used  only  as  a  guic^e  to  optimize 
the  proportion  of  the  dimensions. 

The  expression  (1)  can  be  interpreted  as  if  a  concentration  factor  of  2 
is  applied  to  uniformly  distributed  stress  <jg  in  the  tangential  direction 

while  ignoring  the  radial  stress  a  .  (See  Figure  10^) 

r  ^  • 


Figure  10.  Joint  Design  Configuration  After  First  Tensile  Load. 


In  some  cases,  the  design  of  the  bolted  joint  is  based  upon  not  only 
ultimate  loads  but  also  maximum  deflection.  For  a  joint  designed  accord¬ 
ing  to  the  expression  (1)  (t/a  *  1),  the  hole  elongation  for  glass-epoxy 

joints  is  found  to  be  about  a/2  in  the  instant  previous  to  failure.  For 

other  load  values,  the  deflection  can  be  reasonably  estimated  by  assuming 
a  linear  load-deflection  curve. 
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RESULTS  OF  TESTS  ON  BOLTED  JOINTS:  OPTIMUM  t/a  RATIO 


A  series  of  tests  were  performed  to  determine  the  optinum  t/a  ratio.  The 
bolted  joint  configuration  used  for  this  purpose  is  shown  in  Figure  11. 


(All  Dimensions  in  Inches) 


Figure  11.  Joint  Configuration. 
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The  test  results  are  summarized  in  Table  I.  From  the  values  in  this  table, 
the  ultimate  failure  load  is  plotted  as  a  function  of  t/a  in  Figure  12.  As 
this  figure  indicates,  the  optimum  t/a  ratio  is  approximately  unity. 


TABLE  I 
TENSILE  TESTS 


RADIUS  a 

ULTIMATE  FAILURE 

AVERAGE  OF  ULTIMATE 

SPECIMEN 

INCHES 

RATIO  t/a 

LOAD  103  LB 

FAILURE  LOAD  103  LB 

1 

0.375 

1.67 

38.6 

2 

0.375 

1.67 

39.3 

38.8 

3 

0.375 

1.67 

38.4 

4 

0.437 

1 .29 

43.5 

45.1 

5 

0.437 

1.29 

46.7 

6 

0.500 

1 .00 

54.2 

7 

0.500 

1.00 

55.0 

54.3 

8 

0.500 

1.00 

53.6 

9 

0.563 

0.78 

52.7 

53.1 

10 

0.563 

0.78 

53.4 

11 

0.625 

0.60 

50.8 

12 

0.625 

0.60 

49.7 

50.1 

13 

0.625 

0.60 

49.9 

Figure  12.  Ultimate  Load  as  Function  of  t/a 


MULTIPLE -BOLTED  JOINT 


Many  structures  require  the  use  of  multiple-bolted  joints,  rather  than 
joints  for  only  one  bolt.  With  the  multiple  joint,  the  force  transmission 
occurs  over  a  more  extended  region,  thus  avoiding  the  load  concentration 
on  a  reduced  area  which  occurs  with  the  simple  joint.  However,  in  general, 
the  specific  strength  of  the  simple  joint  is  greater  than  that  of  the 
multiple  joint  when  composite  materials  are  employed. 

The  design  of  multiple,  joints  is  completely  analogous  to  that  of  simple 
joints.  All  the  considerations  given  in  the  first  section  of  this  report 
are  applicable  to  these  joints.  The  total  tensile  or  compressive  load 
must  be  divided  in  equal  parts  between  the  simple  component  joints,  each 
of  which  can  then  be  individually  treated  as  the  simple  joint  of  the 
previous  sections. 

Figure  13  shows  the  multiple  joint,  fabricated  and  tested.  For  the  fabrica¬ 
tion  of  this  joint,  the  joint  of  Figure  14  was  made  using  the  winding  guide 
of  Figure  15  and  the  press-mold  of  Figure  16.  After  molding,  this  joint  was 
cut  into  three  parts  and  bonded,  as  indicated  in  Figure  13,  to  obtain  the 
multiple  joint.  The  simple  tape  layers,  with  fibers  oriented  perpendicular 
to  the  loading  direction,  have  no  principal  structural  function,  but  they 
assist  the  bonding  by  holding  the  parts  together.  Multiple  joints  fabri¬ 
cated  by  this  procedure  are  shown  in  Figure  18. 

The  fabricated  multiple  joints  were  tested  in  tension,  using  the  attach¬ 
ments  shown  in  Figure  19.  Figure  20  shows  the  multiple  joint  mounted  in 
the  test  machine.  Two  specimens  were  tested  in  tension.  The  ultimate 
failure  loads  were  23,800  and  24,800  lb,  respectively.  Figures  21  and  22 
illustrate  the  failure  mode. 

The  failure  load  for  this  type  of  joint  is  estimated  using  the  expression 
(1),  obtaining 


Pult  *  180>000  x  oa  x  °-25  x  6  “  27,000  lb, 

where  a  ^  *  180,000  psi  is  assumed  for  the  glass-epoxy  composite  and  a 

perfect  uniformity  in  the  load  distribution  is  assumed  for  each  bolt.  It 
must  be  noted  that,  for  these  joints,  the  t/a  ratio  is  0.5  instead  of  the 
optimum  t/a  *  1.  These  dimensions  were  adopted  for  construction  reasons. 
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4 


Figure  13.  Multiple  Joint,  Final  Configuration. 
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(All  Dimensions  in  Inches) 


Figure  14.  Multiple  Joint,  First  Stage  of  Fabrication 


Figure  15.  Winding  Guide  for  Multiple  Joint  Fabrication. 


Figure  17.  Multiple  Joints 


QUANTITY:  2  (STEEL) 


PIN  3/16,  LENGTH  2.5,  STEEL  SAE  1050 

QUANTITY :  14 

Figure  19.  Multiple  Joint,  Test  Device. 
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Figure  22.  Tested  Multiple  Joint. 
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IN- PLANE  BOLTED  JOINT 


This  section  considers  a  joint  with  the  bolt  in  the  loading  plane. 

The  configuration  of  this  joint  is  shown  in  Figure  23,  where  only  half 
of  the  fabricated  joint  is  represented.  Figures  24  and  25  are  drawings  of 
the  winding  device  and  the  press-mold,  respectively,  utilized  in  the  fab¬ 
rication  of  the  specimens. 


Two  different  failure  modes  can  be  expected  in  this  joint.  Because  the 
t/a  ratio  is  only  0.333,  a  fiber  failure  of  the  type  shown  in  Figure  6  is 
one  of  the  possible  modes  of  failure.  Another  is  the  transverse  failure 
shown  in  Figure  26,  where  a  debonding  of  the  wound  part  occurs  along  line 
AB.  The  debonding  is  due  to  the  tensile  stresses  acting  perpendicularly 
to  the  loading  direction.  The  stress  analysis  of  the  joint  shows  that 
this  stress  cr^  acting  across  AB  is  only 


a 


n 


_Za. 

43 


where  crn  i-S  the  tangential  stress  at  CD.  Thus,  if  a  reaches  the  value 

180,000  psi,  a  is  about  4,200  psi.  Since  the  transverse  tensile  stress  is 
n 

smaller  than  the  failure  stress  of  the  unidirectional  composite  in  trans¬ 
verse  tension,  this  joint  would  be  expected  to  fail  at  CD  rather  than  across 
AB,  However,  the  failure  occurred  at  AB  in  both  joints  tested.  This  can 
be  explained  if  the  degradation  of  the  composite  tensile  strength  in  the 
fiber  direction  is  considered.  High  tensile  stresses  in  the  fiber  direction 
produce  micro-cracks  in  the  matrix.  These  cracks  result  in  high  stress 
concentrations  in  the  material.  Then,  a  crack  propagation  phenomenon  is 
initiated  along  AB,  and  the  failure  mode  will  be  that  indicated  in  Figure 


26,  even  when  cr  is,  as  an  average  stress  (without  considering  the  stress 
concent  rations  produced  by  the  micromechanical  cracks),  very  small. 


Figure  27  shows  the  joint  mounted  in  the  testing  machine  for  a  tension  test. 
Note  the  attachment  used  in  the  testing.  TWo  joints  were  tested,  both 
failing  along  AB  (Figure  26)  at  loads  of  22,000  and  24,800  lb.,  respectively. 
Figures  28  and  2y  show  a  tested  specimen. 
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(All  Dimensions  in  Inches) 


$3/8 


R  0.3 


n  dc 


n  in  ft  dc 


«Tl 


BUSHING 

(Structural  Aluminum) 


INSERT 

(Laminated  Composite) 


TENSILE 
ELEMENTS 
(Wound  Fiber 
Glass) 


Figure  23.  In-Plane  Bolted  Joint,  General  Configuration. 
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(AIL  Dimensions  in  Inches) 


Figure  24.  In-Plane  Bolted  Joint,  Winding  Device. 


(All  Dimensions  in  Inches) 


MATERIAL:  ALUMINUM 
QUANTITY:  2 


Figure  25.  In-Plane  Bolted  Joint,  Press-Mold. 
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BEFORE 

FAILURE 


AFTER 

FAILURE 


Figure  26.  Failure  Mode. 
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Figure  27.  In-Plane  Bolted  Joint  Mounted  in  the  Testing  Machine. 


Figure  28.  In-Plane  Bolted  Joint  After  Testing. 
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Figure  29.  In-Plane  Bolted  Joint  After  Testing 
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OFF-PIANE  BOLTED  JOINT 


In  this  type  of  joint,  the  axis  of  the  bolt  is  not  on  the  plane  of  loading. 
This  joint  differs  from  the  joints  previously  considered,  in  that  it  in¬ 
cludes  an  adhesive  surface  which  has  a  structural  function  to  perform. 

The  off-plane  joint  is  depicted  in  Figures  30,  31  and  32.  In  the  fab¬ 
rication  of  this  type  of  joint,  the  winding  guide  of  Figure  33  and  the 
press-mold  of  Figure  34  were  used.  The  fabricated  specimens  are  shown  in 
Figure  35. 

Figure  36  shows  the  specimen  mounted  in  the  testing  machine  for  a  tension 
test.  Figure  37  is  the  drawing  of  the  attachment  used  for  the  testing. 

Two  specimens  were  tested,  and  failure  loads  of  2640  and  2915  lb.,  respec¬ 
tively,  were  obtained.  The  failure  mode  is  shown  in  Figures  38  and  39. 

This  failure  mode  is  surprising,  since  shear  failure  along  the  entire  ad¬ 
hesive  layer  (Figure  40)  was  expected.  Instead,  the  failure  occurs  in  a 
region  of  small  tensile  stress  acting  perpendicular  to  the  plane  of  the 
adhesive.  The  failure  shown  in  Figures  38  and  39  does  not  occur  in  the 
adhesive;  rather,  it  is  a  transverse- tensile  failure  in  the  composite. 

This  failure  can  be  interpreted  using  the  arguments  of  the  previous  sec¬ 
tion;  that  is,  it  is  related  to  the  crack  propagation  phenomenon  associated 
with  the  appearance  of  micro-cracks  in  the  matrix. 

From  analysis  of  this  joint  and  the  analysis  stated  in  the  previous  section, 
the  following  conclusion  can  be  made:  in  joint  design  using  composite 
materials,  it  is  necessary  to  eliminate  all  regions  in  which  the  unidirec¬ 
tional  composite  is  under  transverse  tensile  stresses.  Failure  occurs  even 
for  very  small  values  of  these  stresses,  because  of  this  crack  propagation 
phenomenon . 
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26 


3 


Figure  33.  Off-Plane  Joint,  Winding  Guide. 
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(All  Dimensions  in  Inches) 


Figure  34.  Off-Plane  Joint,  Press-Mold. 
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1-0 


Figure  36.  Off-Plane  Joint 


Figure  39.  Off-Plane  Joint  Specimen  After  Testing. 


Adhesive  Layers 


Figure  40.  Adhesive  Layers. 
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ORTHOGONAL  JOINT 


This  type  of  joint,  shown  in  Figure  41,  is  used  to  introduce  small  forces 
into  a  wall.  It  facilitates  the  application  of  handles  or  other  devices. 
The  drawing  of  the  fabricated  joints  is  shown  in  Figure  42,  while  Figure  43 
shows  the  unidirectional  composite  ring  from  which  parts  of  the  joint  are 
obtained.  Figure  44  is  the  winding  device  utilized  to  fabricate  the  ring 
of  Figure  43. 

Wall 


Figure  41.  Orthogonal  Joint. 


Figure  45  shows  three  fabricated  specimens.  Four  specimens  were  tested  in 
the  manner  shown  in  Figure  46.  Figure  47  is  the  drawing  of  the  test  sup¬ 
port  used  in  this  test.  The  failure  loads  were  1845  and  1970  lb.,  respec¬ 
tively,  in  the  first  two  specimens;  while  in  the  other  two,  failure  loads 
of  2975  and  3005  lb.,  respectively,  were  obtained.  The  only  difference 
between  the  two  sets  of  specimens  was  the  adhesive  specification:  for  the 
first  two,  Narmco  227  adhesive  was  used;  for  the  other  two  specimens,  an 
epoxy-polyamide  adhesive  was  used.  The  results  obtained  in  the  testing  of 
the  last  two  can  be  considered  satisfactory.  Figure  48  shows  bonding  fail¬ 
ure  in  the  first  two  specimens.  Figure  49  shows  the  failure  of  a  specimen 
fabricated  with  the  epoxy-polyamide  adhesive.  In  this  case,  the  adhesive 
bond  is  improved  and  the  failure  occurs  in  the  composite  itself,  as  well  as 
in  the  adhesive  layer. 
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7.0 


(AIL  Dimensions  in  Inches) 


Figure  42.  Orthogonal  Joint,  General  Configuration 
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(All  Dimensions  in  Inches) 


1 .25 


Wound  Fiberglass, 
Fiber  Direction 


Parts  to  be  utilized 


1.25 


Figure  43.  Orthogonal  Joint,  Fabrication  of  Parts. 
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MATERIAL:  ALUMINUM  *3/4 


Figure  44.  Orthogonal  Joint,  Winding  Device. 
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Figure  45.  Orthogonal  Joint  Specimens. 


Figure  46.  Orthogonal  Joint  Mounted  in  the  Testing  Machine 
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(All  Dimensions  in  Inches) 


Figure  47.  Orthogonal  Joint  Test  Attachment. 
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Figure  49.  Orthogonal  Joint  With  Epoxy- Polyamide 
Adhesive  After  Testing. 
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LAP  JOINT 


The  lap  joint  is  one  of  the  most  important  types  of  joints,  especially  in 
structures  of  composite  materials,  where  the  use  of  rivets  and  bolts  is  not 

ad v isab le . 

The  correct  design  of  lap  joints  requires  an  understanding  of  the  mechanism 
Lor  load  transfer  from  one  part  to  the  other  through  the  adhesive.  The 
shear  stress  distribution  was  analyzed  in  References  1  through  4  under  the 
assumption  of  linear,  perfectly  elastic  behavior  of  both  adhesive  and 
adherent.  This  study  shows  that  the  shear  stress  concentration  is  very 
high  for  the  joint  dimensions  usually  required  by  the  technical  applications. 
On  the  other  hand,  most  of  the  adhesives  in  common  use  exhibit  a  nonlinear 
stress-strain  relationship,  and  may  sustain  considerable  strain  before  fail¬ 
ure.  The  behavior  of  the  adhesive  can  still  be  considered  perfectly  elastic, 
in  the  sense  that  no  residual  deformation  remains  after  the  unloading,  even 
from  stress  levels  approaching  the  failure  stress.  A  typical  stress-strain 
curve  for  an  adhesive  is  shown  in  Figure  50.  A  convenient  analytical  expres¬ 
sion  for  such  curve  is: 


figure  50.  Typical  Stress-Strain  Curve  for  the  Adhesive, 


Figure  51#  Lap  Joint  General  Nomclature. 
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(2) 


T  -T 


whereT  .  and  k  are  constant.  The  initial  shear  modulus  G  is 


o  \d  y/. 


kT, 


(3) 


Y  =  0 


Let  us  consider  the  lap  joint  shown  in  Figure  51.  Applying  the  equilibrium 
condition  to  the  section  depicted  in  Figure  52, 


'1 


dx 


L 


-T  dx 


.T  dx 


□ 


P8+dP, 


Figure  52.  Infinitesimal  Element, 
the  following  equations  are  obtained: 


dpx 

T  *  9  Px+Pa  *  P  *  constant 

/ 

The  shear  deformation  in  the  adhesive  is 


(4) 


Y  = 


U1  _  u2 


(5) 


where  u^  and  U2are  displacements  of  the  two  adherends. 

Differentiating  y  of  equation  (5)  with  respect  to  x,  the  following  relation 
is  obtained* 


®i  “  ea  *  t 
1  *  a  dx 


(6) 
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Assuming  Hooke's  law  to  be  valid  in  the  adherend, 


Est,  * 


(7) 


equation  (6)  can  be  rephrased 


dV 


(8) 


Ej  Eg  tg  a  dx 
Differentiating  (8)  with  respect  to  x  and  applying  equations  (A),  we  find 


T(v) 


1 


/  i  i 

4  W 


^El+  *,0  (•  k  dxC,| 


ft  ^  +  £ 

la  dx  *  tgEg 


\  —  (—  +  -i-\ 

/  dx  ta^/ 


(9) 


If,  in  particular,  tx  and  tg  are  constant,  this  equation  becomes  the  familiar 
expression  for  the  double  straight-lap  joint: 


t(y)  •  — l 


7~T  +  77  dx 

h 


(10) 


Finally,  introducing  the  definition  (2),  a  nonlinear  differential  equation 
in  \  is  obta ined  : 


£l. ,  (i-«'kY) 

i  2  l  '  ' 


1  1 
+ 


<il) 


dx 

a 

Equation  (11)  is  a  nonlinear  differential  equation  of  the  form 


f(v) 


and  will  be  integrated  numerically. 


ne  boundary  conditions  are 


or,  from  equations  (4)  and  (8), 


(±L\  _ _ £_  (±L)  .  ■  ■  ?■-, 

Vdx/  ta  t  E  *  \dx/  .  t  E 

x  ■  o  3  a  x  ■  L  1  a 


(12) 


Therefore,  the  problem  of  finding  the  shear  stress  distribution  in  a  lap 
joint  with  nonlinear  adhesive  requires  the  solution  of  the  differential 
equation  (11)  with  the  boundary  conditions  (12)* 

The  function  that  solves  the  differential  equation  (11)  has  the  form  shown 

in  Figure  53.  The  values  yf  and  yj  of  the  derivatives  at  the  ends  are  known 
(equations  (12)),  but  the  values  y  and  y  are  not.  This  problem  must  be 

O  v 

solved  using  a  trial-and-error  method. 

Let  us  assume  certain  initial  values  Y^/2  ar)/*  Y ^or  t*ie  ^unction  Y  an<* 

l 

its  derivative  in  x  ■  Applying  the  Runge-Kutta  formulas  to  equations 
(11),  the  strains  and  their  derivatives  are  found  at  the  ends;  i.e., 


yi/l 


\?2 


(13) 


Incrementing  the  values  of  and  by  AY{  /2  and  ^\/2’  resPectively> 

the  effects  of  these  increments  on  the  end  values  can  be  determined: 
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V2  +  A^(/2  ’  yt/2 - Runge-Kutta  ,  y™ ,  y™ ,  Y^° 


Y(l) 

yu? 


'(1).  *  '(1) 


7  /2 


+  Ay 


1/2 


Runge-Kutta 


U)  'CD  (1)  '(1) 

►Yoy  ’  Yoy'  ’  Y(-y '  *  Y-ty ' 


(14) 


With  the  results  of  equations  (13)  and  (14),  the  following  pa  tial  deriva¬ 
tives  "an  be  computed  approximately: 


t 

dy 

To  . 

'(1)  '(1) 

Y  oy  •  Yo 

9Yo 

1— 1 

O 

1 

>- 

•  | 

^/2  " 

^1/2  ’ 

*y't/2 

CM 

O 

1 

SY^  . 

'(1)  '(1) 
Y<y 

'<«  -  v'd) 

.  Y-ty'  yl 

L\n 

by' 

1/2 

</2 

(15) 


On  the  other  hand,  the  increments  of  the  end  values  y^  and  y^  are  related 
to  the  increments  of  the  starting  values  by  the  relations 

3Y'_  By' 


Av“’*vt/2lY«2  +  377^lvi/2 


*>' - 4Vt/2  + 


l  ay 


a  2 


^/2 


Ay4/2 


(16) 


By  taking 


'  ’(1)  t  '(1) 

Y  "  Y  -  Ay0 
o  o 


*  Y£  "  yl 


'(1)  -^(1) 


(17) 


and  combining  with  equations  (16),  we  have 


SY, 


by 


(i) 


by 


AY 


1/2 


1/2  +  ~ 


AY 


1/2 


1/2  “  Yo  '  Yc 


(18) 


-v  ./(D 


SY 


1/2 


V/2  +1y^4y<./2  -  V*  -  Yt 
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Once  the  partial  derivatives  are  computed  using  equations  (15),  the  system 
(18)  cun  be  solved  for  the  increments  Ay^/2  and  /2*  With  these  values, 

the  improved  starting  values  can  be  obtained; 

(2)  (1)  A  '(2)  '(1)  ,  AV' 
yt/2~yl/2  +  L\/2  ’  yl/2  ~  yl  /2  *  L\/2 

(2)  '(2) 

and  the  process  is  repeated.  New  increments  and  ^yl/2  are  adoPted> 

the  partial  derivatives  (15)  are  computed  and  the  system  (18)  is  solved  for 

I 

Ay^/2  and  Ihe  Process  is  considered  complete  when  the  error  in  the 

satisfaction  of  the  boundary  conditions  is  smaller  than  some  value  e.  The 
function  T(x)  is  then  obtained  by  substitution  of  y(x)  into  equation  (2). 

Figure  54  shows  the  stress  diagrams  corresponding  to  a  lap  joint  with  * 

ta«0.0625  in.,  t  =0.004,  £=1.0  in.  The  adhesive  stress-strain  relationship 
is  3 

t(Y)  -  5000  (l-e"36Y)  psi, 

which  corresponds  to  an  initial  modulus 

(G)  =kT=0.18  x  106  psi. 

y=o  £  r 

From  inspection  of  the  curves,  it  can  be  noted  that  the  stress  concentra¬ 
tions  at  the  ends  decrease  when  the  load  is  increased.  This  is  due  to  the 
non-linear  behavior  of  the  adhesive. 

The  ultimate  load  of  the  joint  can  be  evaluated  if  a  design  ultimate  distor¬ 
tion  Y^^  of  the  adhesive  is  specified.  For  example,  for  Y  ^  =  0.04,  the 

ultimate  load  would  be  P  *  1000  lb. 


Figure  55  shows  the  stress  diagrams  corresponding  to  a  joint  with  a  double 
thickness  of  adhesive  (ta  *  0.008  in.).  It  can  be  observed  that  the  stress 
concentrations  are  decreased  below  those  for  the  single- thickness  adhesive. 

Figures  56  and  57  show  the  stress  diagrams  corresponding  to  two  lap-joint 
combinations.  In  both,  one  adherent  is  a  glass-epoxy  composite  while  the 
other  is  either  aluminum  or  steel.  In  these  cases,  the  stress  diagrams  are 

l 

not  symmetric  with  respect  to  x  = 


Four  lap-joint  specimens  were  fabricated  and  tested.  Figure  58  shows  these 
joints.  Figure  59  shows  one  of  them  after  testing.  Table  II  indicates  the 
test  results: 
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SHEAR  STRESS 


SHEAR  STRESS,  psi 


X,  inches 


Figure  55.  Effect  of  Bond-Line  Thickness  t  on  Adhesive 
Shear  Stress  Distributions,  t  -  0.008  inch, 

Ei  *  Eg  =  10  x  10  psi,  t;  =  ta  =  0.0625  inch, 
l  -  1.0  inch,  \  =  5,000  psi,  k  =  36. 
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SHEAR  STRESS,  psl 


SHEAR  STRESS,  psi 


X,  inches 


Figure  57,  Effect  of  Doubler  Modulus  on  Adhesive 
Shear  Stress  Distributions,  Ei  =  30  x  106 
psi  (Steel)  Eg  =  7  x  106  psi,  t*  =  tg  = 
0.0625  inch,  ta  =  0.004  inch,  t  =  1.0  inch, 
»  5,000  psi,  k  =  36. 
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Figure  59.  Failure  in  Straight-Lap  Titanium- 
Titanium  Joint. 
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TABLE  II 

ADHESIVE  LAP-JOINT  TEST  RESULTS 

Specimen 

Failure  Load 
(lb) 

Average  Failure  Stress 
(psi) 

1 

5010 

5010 

2 

5435 

5435 

3 

5780 

5780 

4 

4095 

4095 
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lMlOTOKLAS  TIC  TEST  KESULTS 


In  order  to  evaluate  stress  concentrations  in  bolted  joints,  photoeiastic 
models  were  prepared  and  tested.  Figures  60a  through  60h  show  the  isochro¬ 
ma  tic  lines  and  the  isoclines  (0°,  15°,  30°,  45°,  60°,  75°,  90°)  obtained 
in  a  bolted  joint  under  tension.  It  can  be  observed  from  the  position  of 
the  isoclines  that  the  polar  shear  stress  is  very  small.  This  result 

is  coincident  with  that  obtained  through  finite-element  analysis. 


Figures  61  and  62  show  the  isochromatic  lines  in  a  transparent  isotropic 
model  and  in  the  corresponding  orthogonal  joint  made  from  composite  material. 
In  the  latter,  it  is  possible  to  observe  the  adhesive  bond  lines. 


Figure  60a.  Photoelastic  Model  of  a  Bolted  Joint 
Loaded  in  Tension;  Isochromatics. 


Figure  60b.  Photoelastic  Model  of  a  Bolted  Joint 
Loaded  in  Tension;  0°  Isoclines. 
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Figure  60c.  Photoelastic  Model  of  a  Bolted  Joint 
Loaded  in  Tension;  15'  Isoclines. 


Figure  60d .  Photoelastic  Model  of  a  Bolted  Joint 
Loaded  in  Tension;  30°  Isoclines. 
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Figure  60e.  Photoelastic  Model  of  a  Bolted  Joint 
Loaded  in  Tension;  45°  Isoclines. 


Figure  60f .  Photoelastic  Model  of  a  Bolted  Joint 
Loaded  in  Tension;  60°  Isoclines. 
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Figure  60g.  Photoelastic  Model  of  a  Bolted  Joint 
Loaded  in  Tension;  75°  Isoclines, 


Figure  60h.  Photoelastic  Model  ot  a  Bolted  Joint 
Loaded  in  Tension;  90°  Isoclines. 
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Figure  61.  Orthogonal  Joint  Loaded  in  Tension; 
Transparent  Model. 


Figure  62.  Orthogonal  Joint  Loaded  in  Tension; 
Composite  Material  (Coated). 
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FINITE- ELEMENT  METHOD  FOR  STRESS  ANALYSIS  OF  JOINTS  IN  ANISOTROPIC  MATERIALS 


In  a  previous  study,  a  computer  program  for  the  stress  analysis  of  aniso¬ 
tropic  joints  was  presented.  The  displacement  field  in  a  triangular  element 
was  assumed  linear,  or,  equivalently,  the  stress  field  was  assumed  constant. 

In  order  to  improve  the  accuracy  of  the  results,  a  new  program  was  developed, 
based  on  the  assumption  of  quadratic  displacements  functions.  These 
functions  permit  linear  variation  of  stresses  along  each  element.  In  this 
way  it  is  possible  to  calculate  the  stress  at  a  node  of  the  triangular  net 
simply  by  averaging  the  stresses  corresponding  to  each  triangle  whose  vertex 
is  at  the  node.  This  operation  is  performed  automatically  by  the  program. 

The  accuracy  of  the  stresses  obtained  by  this  averaging  process  is  highly 
satisfactory,  even  on  the  boundaries.  Figure  63  represents  the  basic 
triangular  element  referred  to  the  local  axis  £,17. 


Figure  63.  Basic  Triangular  Element. 


The  coordinates  of  the  middle  points  a»  0>  Y  are  obtained  from 


?b  +  5c 

T|  _ 

v\ 

2 

>  U  * 

a 

2 

5c+5a 

\  +  \ 
c  a 

CM 

’  V 

2 

<a+Sb 

\*\ 

2 

■  V 

2 

(19) 


Thus, only  the  vertex  coordinates  need  be  introduced  as  inputs. 

The  displacement  field  in  the  triangular  element  is  assumed  to  be  given 
by  the  quadratic  functions 


u«  a1  +  a^5  +  33!]  +  a4§^  +  ag^T]  +  a^  Tp 
v=  a,  +  %  5+  *9!)+  a105»+  axl*|)|.  alalf 
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where  through  are  constants. 

Let's  introduce  the  vectors  and  jaj  by  the  expressions 

{6}  ■  {“.  ub  uc  ua  “6  UY  va  vb  vc  va  Y  \  } 

■  fl  a2  a3  a4  a5  a6  a7  a8  a9  a10  all  al2} 


(21) 


Through  use  of  equation  (20),  a  relationship  between  {  6|  and  |a|  can  be 
written: 


{  6  |  ■  [A]  |  a  | 


(22) 


with 


[A] 


[A*]  [0] 
[0]  (A*] 


;  [A*] 


0 

5b 

5, 


l 
l 
l 
l 
1 

Li  s. 


a 


0 

\ 


a 


5 


a  a 


0  1 

2 


5b  5b\  \ 
£  §CT1C  lfc 


£  SA  t 


1f 


5  7|  T f 

Y  Y  Y  -1 


(23) 


From  equations  (22)  and  (23),  the  coefficients  a^  are  determined  immediately: 


{a}  -  [A]"1  jfij  ,  [A]"1  = 


The  strains  in  the  element  are: 


Su 


9v 


[A*]”1  [0] 
[0]  [A*] 


5u  dv 


e  ■  —  »  e  *  rr-  >  y  *  rr:  +  rr 
x  5?  y  ST]  xy  91]  55 

or,  taking  into  account  equations  (20), 


(24) 


V  a2  +  2  a4§  +  a5T) 

V  a9  +  2  ^  +  ^ 


xy 


S  +  ar.  5  +  2  afiTl  +  aft+  2  a^S  +  a^T) 


11 


(25) 
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I*Vr  via  anisotropic  material,  the  stress-strain  relationship  is  expressed  as 


a  =  C  e  +  C10e  +  C1 _Y 
x  11  x  12  y  13  xy 


O'  “  +  C.j,, £  *+• 

y  12  x  22  y  23  xy 


(26) 


T  —  C  ~  ~  ^  4  C ^  4  C0«S 

xy  33  xy  13  x  23  y 

The  strain  energy  for  the  element,  using  equations  (25)  and  (26),  becomes: 

W  =  '  J'J'1  {cn(a2  +  4j4^"  +  a5^  +  ^a2a4^  +  2a2a5T'  +  ^®4a5^) 

+  2Cl2(a2a9  +  Vll*  +  2a2a!2T1  +  2V9§  +  Vll*2 

+  4a4ai2CT1  +  a5a9T1  +  a5all?T1  +  2a5ai2T|2) 

+  2C13  (a 2a3  +  Vs*  +  2a235T1  +  a2a8  +  2VlO?  +  Vll*1 

+  2a3a4§  4  2a La^2  4-  4*^1)  4-a^a^  +  ^ain§2  (27) 


4  5 


4  6  '48s 


4  10 


+  2a4au^  +  a3a5^  +  a5^  ■*"  2a5a6^'^  +  a5a8^  +  2a5ai0^ 


+  a5ail 


Tl2) 


+  2C22(a2  4  ax2  f2  +  4  2a9an§  4  4aga12T]  4-  4a11a12?Tl) 

+  2C23  (a  3a9  +  a5a9^  +  2a6a9T1  +  a8a9  +  2a9al0?  +  a9allT1  +a3allC 
4  a5a11?  +  2aea 1 1^1^1  +  a8all^  +  2a10allF  +  all  ^ 


4  2a3ai2T]  4  la^a^l^  +  4a5a121l  +  2aga12T|  4  4a10a12?1l 

4  2ana12Tl2  ) 
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+  c33(*3  +  *5  ?Z  +  ‘^2  *  a8  +  4a10f2  +  al^  +  2a3a5? 

+  +  2a3ag  +  ^a3axo^  +  2a3ail^  +  4a5a6P^ 

+  2a5a8^  +  4a5aiOc2  +  2a5ail?T]  +  4a6a8T1  +  8a6alO^T| 

+  4a6all'  +  4a8a10^  +  2a8all^  +  4alOall^)  } 

The  variation  of  the  thickness  t  is  assumed  to  be  linear  in  the  element. 
Thus , 

t  ■  Tj§  +  T2T1  +  T0  ,  (28) 


where  the  constants  T^, 

t  at  the  three  vertices: 
c 


and  T 

0 


can  be  found  from  the  thicknesses 


t  ,tw 
a*  b 


9 


(29) 


The  equations  of  the  sides  of  the 
triangular  element  depicted  in 
Figure  64  can  be  written  as: 


,  He 

R2j:Tl-m2?  ;  m2-  — 

c 

R3):T|  -  m35  +  b3;  m3  ■ 


V\ 


.  bo  ■  T)*.  -  mo§. 


.  ,  “3  '.b  -3»b 


B 


Figure  64.  Triangular  Element. 
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Let's  define  now  the  following  integrals,  which  will  be  used  to  express  the 
strain  energy: 


64 


With  equations  (28),  (29),  (30)  and  the  integrals  y^  through  y^,  the  strain 
energy  (27)  can  be  expressed  as: 

w  -  h  Ja(  T  [B]  |  a  |  ,  (31) 

where  the  matrix  [B]  is  the  one  given  on  the  next  page.  To  simplify  the 
writing  of  [B],  the  following  nomenclature  was  used: 


T,1  ■  T1Y1  *  T2v2  +  Vo 

Ty3  '  T1Y3  +  T2Y5  +  Vl 

T  -  Ty  +  Ty  +  Ty 

Y5  1  5  2  4  0  2 

V  *  Vs  +  T2Y7  -  T0Y3 


(32) 


V  '  V?  +  Vs +  V5 

V  -  Vs  +  T2Y9  +  T0v4 

Introducing  equations  (24)  into  equation  (31),  the  strain  energy  may  be 
expressed  as  a  function  of  {6}: 
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*  * 


4 


C  T 

11  yl 

C13Tyl 

2CllTy3 

C33TY1 

2C13TV3 

MW 
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2C11TY7  +  2C13TY6 


C11TY8  +  2C13^7  +  C33TY6 


* 


13  +  C131\5 

2C12TV5 

f  5  +  C2JTV3 

2C23TV5 

rV6+  2C13Tv7 

4C12Ty7 

+  C33)Tv7  +  C13Ty8  +  C23Ty6 

2C12Ty8  +  2C23Ty7 

rY8  +  2C23TV7 

4C23Ty8 

^5  +  C23Ty3 

2C23TY5 

t3  +  C23TY5 

2C22TY5 

ry7  +  2C23TY6 

4C23Ty7 

7 6  +  C33Ty8  +  2C23Ty7 

2C22Ty7  +  2C23Ty8 

4C22TY8 

W  -  h  (5  }T  ([A]'1)T[B]  [a]_1{6  } 

The  nodal  force  vector  fP)  can  then  be  found  by  using  the  second  theorem 
of  Castigliano: 


fp}  -  IqV} =  [Ar^s}, 

and  the  12x12  stiffness  matrix  [K]  of  the  triangular  element  is: 

[K]  =  ( [A]_1  )T[B]  [A]" 1  (33) 


For  certain  boundary  conditions,  it  is  convenient  to  refer  the  stiffness 
matrix  to  rotated  nodal  directions.  Thus,  the  new  stiffness  matrix  [K1] 
is  given  by: 

[K- ]  =  [T]T[K]  [T]  (34) 
where  the  transformation  matrix  [T]  is  expressed  as: 


[T]  - 
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0 
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0 
0 
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0 
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Binra. 
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-  •  trip 

0 

0 

0 

0 

0 

c°«»> 


The  angles  cpa  through  cp  are  the  rotations  at  the  nodes  a  through  Y, 

a  y 

respectively. 


By  adequately  combining  the  matrix  [K']  of  the  elements,  a  total  stiffness 
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matrix  [K*]  for  the  whole  structure  is  obtained.  Reference  5  provides 

details  of  the  procedure  to  find  [K*] .  The  nodal  equilibruim  equation  for 
the  whole  structure  can  then  be  written  as: 


{P*|  -  [K*]  {6*} 


(36) 


It  is  convenient  to  partition  this  matrix  equation  in  the  following  form: 


where  the  symbols  have  the  meaning: 


(37) 


vector  of  unknown  displacements 


vector  of  known  displacements 


vector  of  known  forces 


vector  of  unknown  forces 


From  equation  (37), 


KJ  W-  Kl  -  KJK}. 

{•k 

6  i 

Then,  also  from  (37),  we  obtain: 

It  is  convenient  to  compute  this  vector  |  P£  |  f°r  checking  purposes, 
although  it  is  not  really  necessary  for  the  computation  of  the  stresses. 

Once  all  nodal  displacements  are  known,  the  vector  {a}  for  each  element 
can  be  obtained: 


70 


| a}  -  [Af^f  -  CAj”1  [X] J6 ‘  }  , 

where  |6*  [are  the  displacements  referred  to  the  rotational  directions. 
Finally,  the  stresses  are  given  by: 


ax  "  Cll(a2  +  2a45  +  a5?l)  +  C12(a9  +  all?  +  2al2T1) 
+  C13(a3  +  a55  +  2a  ^T]  +  ag  +  2a1Qf  +  a lJT\) 

Oy  m  ^^2^2  *  2a4^  +  a^)  +  ^22^a9  +  all^  ***  2al2^^ 
+  C23(a3  +  a5f  +  2a&Tl  +  ag  +  2a  1Q?  +  a^) 

T  xy  “  C13(a2  +  2a4^  +  a5T1)  +  C23(a9  +  aU^  +  2a12T|) 
+  C33(a3  +  a5S  +  2a6Tl  +  ag  +  2a1Q?  +  auT|) 


(38) 


The  value  of  the  stresses  at  each  nodal  point  of  the  triangular  net  is 
obtained  by  averaging  the  values  of  the  stresses  at  this  node  from  each  one 
of  the  elements  containing  the  node. 


The  material  constants  C 


ij 


are  obtained  by  assuming  that  the  material  is 


anisotropic  at  each  triangle  and  is  formed  by  fibers  oriented  in  one 
direction  in  the  matrix.  The  axes  x',  y'  in  Figure  65  coincide  with  the 
fiber's  direction,  while  x,  y  are  the  general  axes  of  reference,  parallel 
to  5,  7). 


y 
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The  elastic  constants  with  respect  to  the  x',  y'  axes  are  given  by; 


E  E, 

m  t 


E  ,  -  E  V  +  E,V,  ,  E  ,  ■ 
x  mm  f  f  y  E  V,  +  E  / 

J  m  f  f  m 


G  G, 
m  f 


Vx'y'  "  VmVm  +  VfVf  ’  Gx'y'  “  G  V.  +  G.V 
J  J  m  f  f  m 


where  V  ,  Vr  are  the  volumetric  fractions  of  the  matrix  and  fibers,  E  , 
m  f  m* 

G  ,  v  are  the  elastic  constants  of  the  matrix  material;  and  E£,  Gr,  are 
mm  f  r  f 

the  elastic  constants  of  the  fiber  material.  From  E  • ,  E  f,  Vv  r  i  and 

X  y  x  y 

G  ,  ,  ,the  following  new  set  of  constants  is  found: 
x  y 

.  i  _1_  ■  .  Vx  'y 1  ,  1 _  ,  '  1 _ 

&11  "  E  ,  *  12  *  "  E  ,  ’  D22  E  ,  ’  33  "  G  ,  , 

x  x  y  x  y 


The  elements  of  the  symmetric  matrix  [ C  *  ]  referred  to  the  axes  of  orthotropy 
x ' ,  y '  are: 


b22 

r 1  ^  r  i 

C11  A  ’  C12 


b12  bll  1 

_i£  o'  _  _ii  o'  _  _i_  o'  .o'  _  n 

A  ’22  A  ’  33  b^3  ’  C12  L23 


L  “  bUb22  “  b12 

The  matrix  [G]i  referred  to  the  axis  x  and  y,  is  obtained  from  [ C  1  ]  by 
performing  the  operation 


[cj  -  [Tjfhc'ny, 


with 


m  2  2 

2  2 

£  m  2lm 

K*  m  m 

2  2 

m  t  -2tm 

[T2l  - 

2  -  2  . 

m  L  -tm 

-f,m  K/  m  Kf  —  m 

m  _ 

-2t,m  2tm  •  m ^ 

L  ■  cosg 
m  -  sing 
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It  must  be  noted  that,  in  general,  and  are  not  zero  because  x  and 
y  are  not  necessarily  axes  of  orthotropy. 

The  computer  program  is  contained  in  the  appendix  to  this  report. 
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APPENDIX 


COMPUTER  PROGRAM  FOR  FINITE- ELEMENT 
ANALYSIS  OF  JOINTS  IN  ANISOTROPIC  MATERIALS 


program  sixpt 

FD (  I  )  -  FORCE  OR  DISPLACEMENT  AT  NODE  I*  WHICHEVER  IS  KNOWN 

ANISOTROPIC 

NO(I.J)  -  NODE  TYPE  j  AT  TRIANGLE  I 
r  Nn (  I  ,n  -  A  N<TF 
C  NO (1*2)  -  gamma  NODF 
C  NO(  I  *3)  -  R  NODE 
C  NO(  1*4)  -  ALPHA  NODF 
C  NO<  1 .5  )  -  C  NODE 
C  NO(  1*61  -  BETA  NODE 
C  ETA (  I)  -  ETA  AT  NODE  I ) 

C  XI (  I  )  -  X  I  AT  NODE  I 
C 

C  NT  -  NUMBER  OF  TRIANGLFS 

C  NV  NUMBER  OF  VERTICES 

C  NMP  -  NUMBER  OF  MID-POINTS 

C  NUD  -  NUMBER  OF  UNKNOWN  DISPLACEMENTS 
C  NKD  -  NUMBER  OF  KNOWN  DISPLACEMENTS 
C  NUF  -  NUMBER  OE  UNKNOWN  FORCES 

C  NKF  -  NUMBER  OE  KNOWN  FORCES 

C  NN  TOTAL  NUMBER  OF  NODES 

C  NKFD  -  TOTAL  KNOWN  FORCFS  AND  DISPLACEMENTS  <«  2*NN) 

C 

C  IDIRUtl)  NUMBER  OF  X-D I SPL ACEMENT  AT  NODE  I 

C  IDI R ( I .2)  NUMBER  OF  Y-D I SPLACEMENT  AT  NODE  I 

C 
C 

C  NODF(I.J)  -  J-TH  NODE  TOUCHING  ITH  DIRECTION 

r 

c 

C  NOT (  I  *  J)  -  J-TH  TR I  ANGLE  TOUCHING  NODF  I  (UP  TO  J«8) 

C  NOT (1*9)  -  NUMBER  OF  TRIANGLFS  TOUCHING  I-TH  NODE 
C 

C  NOT(I.J)  INCREASES  MONOTONICALLY  WITH  J 
C 

C  A  MAXIMUM  OF  80  TRIANGLES 
C  250  NODES 

0  500  DISPLACEMENTS.  KNOWN  AND  UNKNOWN 

C  60  ROWS  IN  A  PARTITIONED  S'JBMATRIX 

C  60  KNOWN  DISPLACEMENTS 

C  8  TRIANGLES  COMMON  TO  A  NODE 

C 

C  IF  TWO  VFRTICES  IN  THE  SAME  TRIANGLE  HAVE  THt  SAME  X  COORDINATE 

C  THEN  THEY  MUST  BE  ORDERED  TO  BE  A  AND  B 

C 

COMMON  /l!  NO (80.6).  FDI500).  NODE(S00.2).  XI250).  Y ( ?S0 * . TH  (  250 ) , 
]  N.  KREM.  M,  5(60.60.5).  IDIR(250.2).  ALF(250).  T( 17.12). 

2  AMX<12.12)«  BMXI12.12).  SPI60.20).  EF(BO),  FR(80),  RNU ( 80 ) « 

3  FNUI80).  NOT ( 25 J . 9 ) .  C  I  J  ( 6 .80 ) . KMX  (  1 2 . 1 2 . 80  I .  AS(6.6.80). 

4  G ( 6 0 ) .  NKD.  NUD.  NN.  NT  .  V F  (  80 ) .Bt TA ( 80 ) 

DC  A  I  IfUX 

DIMENSION  DUMU3) 

DEG  *  .0174533 

r 
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C  ZERO  OUT  MATRIX  T  AND  MATRIX  R 
DO  U>  I  •  1.12 
DO  1C  J  •  1.12 
BMX(I.J)  « 

10  T(  I.J)  *  0. 

C 

C  ZERO  OUT  APPROPRIATE  PARTS  OF  12X1?  A-MATRIX 

C 

DO  20  I  *  1.6 
DO  ?('  J  =  1.6 
AMXII+6.J)  *  0. 

20  AMXII.J+6)  «  0. 

READ  1000,  ( DUM (  I  1,1*1,13) 

PRINT  1001.  (  Dl/M  (  11,1*1,13) 

READ  1002.  NT.  NV.  NMP ,  NMD.  NKD 
PRINT  1003.  NT.  NV,  NMP,  NMD.  NKD 
NN  »  NV  *  NMP 
NKF  »  NUD 
NUF  *  NKD 
NKFD  »  NKF  «•  NKD 
C 

C  ZFRO  OUT  ALL  KNOWNS 
C 

DO  30  1*1, NKFD 
FD(I)  *  0. 

30  CONTINUE 
C 

C  READ  TRIANGLE  INFORMATION 

C 

00  40  K  ■  1 .NT 

PFAD1007. I.(NO« I .JI . J-1.6)  »FF(  I),FNU( I ) »FR < I ) .RNU (  I  I  , VF (  I ) .BETA ( I ) 
40  CONTINUE 
C 

C  READ  NON-ZERO  KNOWNS  INTO  FD  AND  PRINT  THEM  OUT 
C 

PRINT  1110 
DO  50  1*1. NKFD 
READ  1115.  J.TEMP 
IF  (J.EO.O)  GO  TO  60 
FD( J)  *  TFMP 
50  PRINT  1120  , J.FDI J) 

60  CONTINUE 
C 

C  PRINT  OUT  TRIANGLE  DATA 

r 

LINE  *  4 

PRINT  1030,1  J.J-1 .6) 

DO  80  1*1 .NT 

IFILINE.LT. 54)  GO  TO  75 
LINE  «  4 

PRINT  1030.1 J.J*1, 6) 

75  LINE  »  LINE+2 

80  PRINT  1033,1  .  (N0( I. J) . J*1 .6  )  ,EF( I ) .FNUI I ) ,ER( I ) ,RNU(  I)  , 

1  VF ( I ) ,  BETA! I  ) 

C 

C  PRINT  FD 

C 

C 

C  READ  IN  VERTEX  INFORMATION 
C 

DO  100  K *  1 »NV 

READ  1003,  I,  IDIR(I.l).  I D I P ( I .2 ) . ALF ( I ) ,  X(I),  Y(I),  TH(I) 

ALFIII  -  ALFI  I  )  *  DEG 
100  CONTINUE 
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C  READ  IN  MID-POINT  INFORMATION 
C 

DO  135  K  ■  ltNMP 

READ  1  006  t  It  IDIRUtDt  I  DIP  (  1 1?  )  tALF  U  ) 

ALF { I )  ■  ALF  <  I )  #  DEG 

C 

C  FIND  X t Y  FOR  MID  POINT 

C 

DO  105  L  •  ltNT 
DO  105  J  ■  lt6 
I F ( I.FQ.NOIL tJ) )  GO  TO  110 
ins  CONTINUE 

PRINT  1050  ,1 
STOP 

110  NOLI  ■  NO(Ltl) 

1F<J.EQ#6)  GO  TO  115 

NOLJP1  ■  NOC Lt J+l ) 

GO  TO  120 

115  NOLJP1  *  NO ( L  1 1  ) 

120  NOLJM1  *  NO(LtJ-l) 

130  V ( T )  «  ( Y ( NOL JM1 )  +  Y  ( NOL JP 1  )  )  / 2  • 

XU)  «  (XINOLJM1  H-X(NOLJDl)  )/2# 

TH(  I  I  *  ( TH( NOLJP1 *  +TH(N0LJM1  )  )  /  2* 

135  CONTINUE 
C 

C  PRINT  OUT  NODE  INFORMATION 
C 

I  -  1 

140  CONTINUE 
LINE  «  4 
PRINT  1020 
145  CONTINUE 

BLD  *  ALFI I ) /DEG 

PRINT  102 1 •  It  I D IR  (  Itl)tIDIR( 1 1 2 ) t  X ( I  )  *  Y  < I )  t  TH  (  l )  tBLD 
I  *  1  +  1 

I F {  NGT.NN)  GO  TO  160 
LINE  «  LINE  +  2 
IF ( LlNEtGT *53 )  GO  TO  140 
GO  TO  145 
160  CONTINUE 
C 

C  READ  IN  PARTITIONING  INFORMATION 

C 

READ  1 0 16 1  MSIZE  ♦  KREM  ,  N 
M  «  MSIZF  /  N 
I F ( KREM«NE*N )  M  «  M+l 
C 

c 

C  GENERATE  TABLE  OF  TRIANGLES  VS  NODES 
C 

DO  200  I*  1 1 NN 
NOT ( I 1 9 )  «  0. 

?00  CONTINUE 

DO  300  I  *  1 1 N T 
DO  300  J  *  1 1 6 
K  -  NO(liJ) 

NOT ( K  t9 )  ■  N0T(K»O)  +  1 
L  *  NOM  K  t?  t 
NOT(KtL)  *  I 
3on  continue 


c 
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C  GENERATE  table  OF  NODES  VS  DIRECTIONS 


DO  350  1*1. NN 
DO  ISO  J«1  .? 

K  *  IDIR( I.J) 

NODE  I  *c  .2  I  *  J 
35n  NOOFIK .  1  )  *  I 

GFNERATf  SMALL  K-MA TRICES  AND  ALL  A-STAR  MATRICFS 


DO  40o  I  *1  .NT 

BETA!  I  )  »  BETA!  I  )*DEC> 

CALL  SMALL*!  I  ) 

400  CONTINUE 

CALL  SOLVIT 

1000  FORMAT ( 1 3 A6  ) 

1001  FORMAT ( lHl»13Aft////  ) 

1  JO?  FORMAT (513) 

1 003  FORMA’’  (IX.  6X.9HTR JANGLFS,7X.BHVrRT  ICES  •.  5X  .  1  OHM  I  0-POINTS  , 

1  5X  .  1  3HUNKN0WN  DISPS.5X.1  1HKN0WN  I)  I  FP  S/  1  X  3  (  1  2X  .  I  3  )  .  1  5X  .  I  3 . 13X  . 

2  13) 

1005  FORMAT  ('M3.4F6.3  ) 

1006  FORMAT ( 313. Eft. 3) 

1007  F0RMATI7I3.6E9.3) 

101ft  FORMAT (313) 

1020  FORMAT (5H1NODE.3X.5HDIR  1.3X.5HDIR  2 . 9X . 1HX ,9X . 1HY  .4X .9HTH ICKNESS . 
1  5X.5HANGLE/I 

1021  FORMAT  (?X.I3.2(5X.13),4F10.3/I 

1030  F0RMAT(9HlTRIANGLE.ft(?X.3HPT  . 1 1 ) .  1  OX .2HFT ,BX . AHNtl  F.10X.2HER. 

1  BX.AHNU  R.10X.2HVF.8X . AHBET A/ I 

1033  FORMAT (6X. I3.6( 3X. 13  )  .6(L12.5 ) /) 

1050  FORMAT (5H1N0DE. I3.21HNOT  IN  TRIANGLE  INPUT) 

1115  FORMAT  (  I3.E15.fe) 

1118  FORMAT  (  1H1  *  1  >  >  X . 28HNON-2ERO  BOUNDARY  CONDITIONS  //) 

1120  FORMAT ( 7H  KNOWN!  .13, 3H)  *, E15.fi) 

STOP 

END 
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SUBROUTINE  LARGE* < ICYCLE .A.B.C *CS*0*SK12 1 1 1  *  12 . I  3 * 14 . 15 ) 

C 

c  ANISOTROPIC 
C 

C  THIS  SUBROUTINE  IS  CALLED  RY  SUBROUTINE  SOLVIT  TO  GFNERATE  A  ROW  OF 
C  SUBMATRICES  FROM  THE  MATRIX  K-STAR.  THE  SMALL  (12X12)  K-MATRICFS 

C  GENERATED  ON  THE  PREVIOUS  PASS  WILL  BE  SAVED  ON  THIS  PASS*  AND  THF 

C  REOUI RED  SMALL  ^-MATRICES  NOT  ON  HAND  WILL  Rt  GENERATED.  THE  MATRIX 
C  ASTAR-INVERSE  WILL  Be  SAVED  FOR  THF  CALCULATION  OF  STRESSES*  AS  WILL 
C  C 1 1 *C12  *C2?  *C33  FOR  EACH  NEW  TRIANGLE. 

C 

C 

C  A  CORRESPONDS  TO  A<  1*1-1)  IN  THF  PARTITIONED  MATRIX 

C  B  CORRESPONDS  TO  A(It!  )  IN  THE  PARTITIONED  MATRIX 

C  C  CORRESPONDS  TO  A(  1*1  +  1)  IN  THE  PARTITIONED  MATRIX 

C  THEREFORE.  A ( 1  )  IS  THE  TRANSPOSE  OF  C(I-l),  I.CT.l 

C 

r 

C  ON  THE  M-TH  PASS*  <2?  STARtNKD  X  NKO )  IS  RFTURNFH  IN  C. 

C 

COMMON  /Z/  NO (60*6)*  FD(500).  NODF(50O*2).  X(250).  Y ( 250 )  * TH ( 2 5 0  )  * 

1  N*  KREM*  M*  S ( 60  *60  *  5 ) *  IDIRI250.2).  ALFI250).  T(12.12). 

2  AMX  (12*12)*  BMX  (  12  *  1 2  )  *  SP(6C*20>*  EF(80)*  LR  (  80  )  •  RNU(80)* 

3  FNU  (  8u  )  •  NUT  (  25  J  *9  )  «  Cl  J  (  6 .80  )  *KMX  (  12  •  12  *00  )  *  AS(6*6.8D* 

4  G<6<>#  NKD*  NUD*  NN.  NT  * VF ( 80 ) *BE TA ( 00 ) 

REAL  <MX 

DIMENSION  A(  11*12)  .  BUl.Il).  CU1.I3)  .  SK12U1.I4) 

DIMENSION  CSU5.I1)*  0U4.I4) 

C 

DIMENSION  RK12( 360U) 

EQUIVALENCE  (S(  7201 ) .RK12) 

DIMENSION  ND I R ( 6*2 )  *  R(80) 

DATA  ND I R  /  1.  6.  ?t  4.  3*  5.  7.  12*  8*  10.  9*  11/ 

ND IRC  I  *  J I  -  POSITION  IN  GENERIC  TRIANGLF  OF  DIRECTION  J  AT  NODF  T 

( J* 1 *2 )  *  ( 1*1*6) 

IG  *  11 
IABC2  *  4 

I F ( ICYCLF^GT.l  )  GO  TO  ID 
REWIND  l5 
I ARC  1  «  2 
IS  *  o 
JS  *  -12 

GO  TO  22 

10  CONTINUE 

I F ( ICYCLE. GT. 2)  GO  TO  15 
JS  *  0 

K 2  =  KREM 

GO  TO  20 

IS  JS  *  KRFW  +  ( I CYCLF-3  ) *N 

20  IS  *  KREM  +  (  I  CYCLE-2  )*N 

IS  -  BIAS  FOR  ROWS 

JS  -  BIAS  FOR  COLUMNS 

I  ABC  1  =  l 

I  F ( ICYCLE. EO.M)  I  ARC  2  *  5 
22  JSA  =  JS 

JSR  *  JSA  +  12 
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ISC  ■  JSB  ♦  I  i 
IROwl  •  1  ♦  IS 
IROW2  ■  11  ♦  IROwl  -1 
DO  l SO  I  ARC  *  I  ABC  1  t  I ABC2 
GO  TO  (2S*35*4J*41.49)*IABC 
?S  REWTND  11 

RFAO  (  11)  (  (  A<  I  *J)  #1»1  *11  )*J»1*!2> 

GO  TO  ISO 
IS  I  COL  1  »  JS6  ♦  1 

ICOL2  «  ICOLl  ♦  II  -  1 
GO  TO  4S 

40  IF(  ICYCLE.EO.M )  GO  TO  ISO 

ICOLl  •  JSC  ♦  1 

ICOL2  •  ICOLl  ♦  13  -1 

GO  TO  4S 

41  ICOLl  »  NUD  ♦  1 

ICOL2  *  NUD  ♦  NKD 

GO  TO  4S 

49  ICOLl  «  NUD+1 

ICOL2  ■  ICOLl  ♦  NKD  -  1 
IROWl  «  ICOLl 
IROW?  •  I COL 2 
4*  CONTINUE 

DO  100  I  *  I  ROW  1 « I ROW2 
NO  I  ■  NODE ( 1*1 ) 

C 

C  NO  I  -  NODE  WHERE  DIRECTION  I  IS  LOCATED 
C 

DO  100  J« I  COL  1 • I  COL  2 
NOJ  *  NODE ( J  1 1 ) 

r 

C  NOJ  -  NODE  WHERE  DIRECTION  J  IS  LOCAIEl, 
C 

TEMP  *  0. 

NOT  1 9  *  NOT (NO  I .9) 

NOT  J9  *  NOT ( NOJ • 9  ) 

C 

c  I  MAX  -  NO  OF  TRIANGLES  AT  NODE 

c 

DO  64  K  *  1 » NOT  I  9 
NOT  I  *  NO  T ( NO  I t  K  ) 

DO  64  KK*  1 tNOT J9 

C 

C  NOT  I  -  KTH  TRIANGLE  TOUCHING  NODE  NO  I 
C 

IF(NOTI.NF*NOT<NOJ*FK)  )  GO  TO  64 
DO  60  l *  1  *6 

I H  NO  I • EQ  *NO ( NOT  I  t  L  )  )  IK  »  L 
I F  (  NOJ* EQ#NO  ( NOT  I  *1.  M  JK  *  L 
60  CONTINUE 

IM  *  NODE ( 1*2) 

C 

C  IM  -  D I  RFC  T I  ON  1  OR  2*  NODE  NOI 

62  NDIRI  *  NHIP (  iKt  IM) 

JM  -  NODF ( J  *  2 ) 

C 

C  JM  -  DIRECTION  1  OR  2t  NODE  NOJ 
C 

NOIRJ  *  NDIRI JK  * JM ) 

TEMP  *  TEmp  +  kmx (  NDIRI  •NDIRJtNOTI ) 

64  CONTINUF 
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6?  GO  TO  (70»70»80*82*85)*I ABC 
70  B ( I - IS* J-JSB  )  »  TEMP 
GO  TO  90 
80  CONTINUE 

C< I-IS.J-JSCI*  TEMP 
CS  (J-JSCtl-IS)  «  TFMP 

CS  IS  THE  TRANSPOSE  OF  C.  SAVED  FOR  USE  AS  A  IN  THE  FOLLOWING  CYCLF 
GO  TO  90 

82  SK12I  I-IS*J-NUD>  •  TEMP 
GO  TO  90 

85  OI I-NUD*J-NUD)  »  TEMP 
90  CONTINUE 

K22-STAR  ( NKD  X  NKD )  LEFT  IN  ISA  ■  0)  AFTFR  M-TH  PASS 
100  CONTINUE 

I F I IABC.NE.3 )  GO  TO  110 
REWIND  11 

WRITE  (11)  ( (CC(J.I).J«ltI5).I-l.Il) 

110  CONTINUE 

I F ( IABC.NE.4 )  GO  TO  150 

CALCULATE  INDEPENDENT  VECTOR  SEGMENT  G 

DO  130  L«1 *IG 
GIL )  »  FDIL-MS) 

130  CONTINUE 

CALL  MX.VULT  (  S  (  1  *  1 » 3 )  .FDINUD+l )  .R.U  .NKD.l ) 

DO  140  L*1 • IG 
G(L)  ■  GIL)  -  RIL) 

140  CONTINUE 
150  CONTINUE 

K21DIM  «  I 1*NKD 

WRITE  111')  IRK12IJJJ)  *  JJ  J*  1  ♦<?  10 1 M ) 


RETURN 

END 
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SUBROUT  INF  SOLVIT 
C 

C  ANISOTROPIC 

C 

C 

C  THIS  ROUT  INF  SOLVES  SU«G  .  WHERE  S  IS  A  TRI-DIAGONAL  MATRIX  IN 
C  SUBMAT R I CF S «  WITH  ELEMENTS  OF  ORDER  N. 

C  S  IS  KNOWN 

C  SP  IS  A  VFCTOR  OF  DIMENSION  INXV)  WHERE  M  IS  THE  NUMRER  OF  DIVISIONS  OF  S 
C  C  IS  WRITTEN  ONTO  TAPE  AFTER  DERIVATION  ON  THE  FORWARD  PASS. 

C  AND  READ  BACK  IN  ON  THt  GACKSWEFP 
C  S(l.l.l)  INITIALLY  CONTAINS  SII.I-1) 

C  S(l.lt?)  INITIALLY  CONTAINS  SII.I  ) 

C  Sll.l. 4)  INITIALLY  CONTAINS  SII.I+1) 

C  SP  CORRESPONDS  TO  P  IN  THE  WRITEUP  BY  GATEWOOD  ON  THE  FORWARD  PASS. 

C  ON  THE  BACKSWEEP.  IT  CORRESPONDS  TO  U. 

COMMON  /Z/  NOIBO.A).  FDI500).  NODFI50O.2).  X<?50),  Y I  250 ) . TH I ?50 J , 

1  N.  KREM.  M,  SI60.60.5).  IDIRI2D0.2).  ALFI250).  TI12.12). 

2  AMXI12.12).  bMX (12.12).  SP ( 60 .20 ) .  EF(80).  ER(80).  RNU(SO). 

3  FNU  ISO).  NOT  I25.I.9)  <  C  I J  (  6 .80  )  .KMX  (  1 2 . 1 2 .80 )  .  ASI6.6.80). 

4  G 1 6  C 1 )  .  NKD.  NUD.  NN.  NT  . VF (SC ) .BE TA  (  80 ) 

REAL  KMX 

EQUIVALENCE  (S< 144011.0 
DIMENSION  Cl  3600) 

DATA  XLIMIT  /l.E-8/ 

REWIND  9 
N2  «  2*N 
KREM2  «  2*KREM 
DO  30  I  CYCLE  »1.M 
I F I  ICYCLE-2 )  1.2.3 

1  K 1  «  KRF* 

K2  »  KRFM 
K3  »  KREM2 
GO  TO  4 

2  K 1  *  N 
K2  ■  KREM 
K 3  «  N? 

GO  TO  4 

3  K1  «N 
K2  »N 
K 3  «N2 

4  CONTINUE 
K4  «  N 

I  F (  ICYCLE.FO.M  )  K<i  *  NKD 

CALL  LARGEKI I  CYCLE . S . S ( 1 . 1 . 2 ) *S ( 1 . 1 .4 ) . S ( 1  ♦  1 .3 ) »S ( 1 ♦ 1 .4 )  , S  ( 1 . 1 .3 ) . 

1  K 1 .K2.K4.NKD.N) 

I F I ICYCLF.EO.l  )  GO  TO  n 

CALL  MXMULTIS(l.l.l)  .  Sll.l. 5)  .  Sll.1.3)  .K1.K2.K1) 

Sll.1.5)  CONTAINS  C  FROM  LAST  CYCLE 

CALL  MXSUBISI1 .1  .2 )  .  Sll.1.3)  .  Sll.l. 2)  .K1.K1) 

10  CONTINUE 

R  (  I  .  I)  NOW  IN  S I  1  .  1 . 2 ) 

CALL  INVFRTISI 1.1.2)  .  K1  K3.  XLIMIT  .  FLAG) 

IF  IFLAG.NF.-. )  GO  TO  600 

INVERSE  OF  BII.il  NOW  IN  Sll.1.2) 

IF  I ICYCLF.FO. 1 )  GO  TO  ?' 
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CALL  mxMULTISI  1.1.1  )  .SP I  1 . ICYCLE-])  ,  Sll.1.1)  ,  N, 
CALL  MXSUB  IG  *  Sll.1.1).  G.  N.  1) 

20  CONTINUE 

CALL  MXMULT I  SI  1 .1 .2 )  .  G  .SP(1.  ICYCLE  )  .  XI  .  Xl 
IF  I ICYCLE. GE.M)  GO  TO  15 

CALL  MXMULT  (Sfl.1.2)  .  Sll.1.4)  .  ,511.1.5).  XI  ,  XI 

NSQ  *  X1*N 

WRITF  (9)  <C( JJJ) .JJJ«1 .NSO) 

10  CONTINUE 
15  CONTINUE 

NOW  IN  RACXSWEFP.  SOLVING  FOR  U 

DO  60  I  *  2.M 

JCYCLE  =  M-I  +  l 

IF I JCYf LE.GT • 1 )  GO  TO  16 

XI  »  XRFM 

GO  TO  37 

36  XI  *  N 

37  CONTINUE 
NSO  ■  X1*N 

IF  ( JCYCLE. EO.M-1)  GO  TO  41 
BACXSPACF  9 
41  CONTINUF 
BACXSPACE  9 

READ  (9)  (C( JJJ) . JJJ-I .NSO) 

U ( M I «SPI  M  )  .  CONSIDER  FIRST  (M-l)TH  CYCLE 

CALL  MXMULT I  SI  1 . 1 .5  )  .SP ( 1 . JCYCLE+1  )  .5(1 . 1 . 1 ) .X 1 «N . 
CALL  MXS'JBISPU.  JCYCLE).  S(l.l.l)  -SP  ( 1  .JCYCLE  )  .  Xl 
60  CONTINUE 

U(NS.l)  NOW  STORED  IN  SPIN. I)  .  I«1.M 

CALL  ENDIT 
RETURN 

500  CONTINUE 

PRINT  1000  ,  ICYCLE 
STOP 

loop  FORMAT  (11H1C0ULD  NOT  INVERT  MATRIX  IN  ROW. I?) 

END 


X?  .  1  ) 

.  1  ) 

.  N) 


1  ) 

.  1) 
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SUB  ROUT  I  ►iF  FND1T 

This  5UPROUT I NF  CALCULATES  the  stresses  sic.max.  SIGMAY,  and  tauxy 
AT  EACH  NODF.  FOR  EACH  TRIANGLE.  THFN  FOR  EACH  NODF.  AN  AVERAGE 
VALUF  is  DERIVED  BY  SUMMING  the  values  arrived  at  IN  EACH 
CONTIGUOUS  TRIANGLE.  AND  DIVIDING  BY  THE  NUMBER  OF  TRIANGLES9 


C  SHIFT  THE  DELS  TO  REMOVE  THE  GAPS 


EQUIVALENCE  (S.SMALLA)  ,  <S< 3601 )>PSTR)  .  (SI10801I.Q) 

DIMENSION  SMALL  A 1 12  •  80)  .  DXII?)  .  SXI3.250) 

DIMENSION  PSTRI 3.250)  .  01 3600) 

EQUIVALENCE  ( S( 7201 ) .SX)  .  ( S ( 8001 ) .DEL ) 

COMMON  III  N0I80.6).  EDI  500  I  *  NODFI500.2).  X(250).  Y ( 250 ) . TH ( 250 1 . 

1  N.  KREM.  M,  SI60.60.5).  IDIRI250.2).  ALFI250).  TI12.12). 

2  AMXU2.12).  BMXt  12.12).  SPI60.20).  EFI80).  ERI80).  RNUI80). 

3  FNU ( DC ) .  NOT  1250.9).  C I J ( 6 .80  )  .KMX ( 12 . 12 .80 ) .  ASI6.6.80). 

4  Glfi”).  NKD.  NUD.  NN.  NT  . VF ( 80 ) .BETA ( 80 ) 

PFA|.  <mx 

PJMFNSION  DFLI5O0)  ,  PS2(6(  )  .  NTRTI6)  .OXP  ( 12  )  .NOD  ( 6 1 

DATA  NTRI/1.3.5,4,6.2/ 

MSI ZE  *  NKD+NUD 
PI  *  3.1415927 
DO  10  J*  1  .KREM 
DEL  I J  )  *  SP(J.l) 

17  CONTINUE 

KLOC  ■  KRFM  -N 
no  20  1*2. m 
. .  KLOC  •  KLOC  ♦  N 
DO  20  J*1,N 
DEL (KLOC+J)  *  SP( J.I ) 

20  CONTINUE 

DO  22  1*1. NKD 
22  DEL C NUD* I )  «  FDINUD+II 
“PINT  10)0 
NDEL  »  ( NUD+NKD )  /  7 
JCNT  *  0 
DO  70  J* 1 .NOEL 
JCNT  *  JCNT  ♦  1 
IFI JCNT.LE.17)  GO  TO  25 
PRINT  101 : 

JC»'T  *  1 

25  JF I  P  «  7*  f  J—  1 )  1 

JL AST  *  JFIR  ♦  6 
PRINT  1011  .  (K.K*JFIR,JLAST) 

PRINT  1012.  ( DEL ( K ) .K* JF IR . JLAST ) 

30  CONTINUE 

LOCI  *  7*NDEL*1 

L0C2  *  NUD  ♦  NKD 

IFI  LOCI. GT.L0C2)  GO  TO  4.7 

PRINT  1011 .(K.K*L0C1 .L0C2) 

PRINT  1012.  (DFLIK) .K*LOC1.LOC2) 

40  CONTINUE 

r 

C  CALCULATE  STRESSES 
C 

DO  100  1*1. NT  ' 


DO  SO  J*l,8 
NTRIJ  *  NTRKJ) 

NODI J)  *  NO  I  I .NTRIJ) 
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NOE  «  NOO(J) 

IDIR1  ■  1 0 1 R  f NOE  *  1 ) 

IDIR2  ■  IDIR (NOE  *2  ) 

DXP(J)  «  DELI IDIR1 ) 

60  DXPIJ+6)  ■  DEL* IDIR2) 

r 

C  CALClJLATF  OrLTA  BY  ROTATING  DXP 
C 

DO  80  J-1.6 
NOE  *  NOD(J) 

CALP  ■  COS <  ALF I NOE  )  ) 

SALP  »  SI N ( ALF ( NOE )  ) 

DXIJ1  ■  CALP*DXP( J)-SALP*DXP< J+61 
80  DXIJ+6)  *  SALP*DXP(J)  +  CALP*DXP( J+6 ) 

CALL  MXMULTI AS< 1 .1 . I  1 .DXC1  1  .SMALLAU . 11,6.6.1) 

CALL  MXMULT  (  AS ( 1  , 1 . I ) .0X17)  »SMALLA(7.I) .6.6.1 ) 

100  CONTINUE 
C 

C  HAVE  VECTOR  SMALL  A  FOR  EACH  TRIANGLE 

r 

C  FIND  SIGMAX.SIGMAY.TAUXY  FOR  EACH  NODF  »  THEN  AVERAGE  AMONG  TR I ANGLFS 

r 

DO  200  I  »  l.NN 
NOT 9  ■  NOT! I .9) 

DO  110  J*1  »3 
SX(J.I)  ■  ■>. 
lio  CONTINUE 

DO  120  J*1 »N0T9 
L  »  NOT! I .J) 

NOL 1  »  NO ( L . 1 1 
XI  ■  XII)  -  X ( NOLI  1 
ETA  •  Y  ( I  1  -  Y(NOLl) 

Cll  ■  CIJI1.L) 

Cl?  *  CIJI2.L) 

C13  ■  CIJI3.L) 

C22  *  CIJ(A.L) 

C23  =  CIJIS.LI 
C33  ■  CTJI6.L) 

A2FTC  «  SMALLA(2.LH-2.*SMALLA(4,LI*XI  +SMALL A ( S  »L  1 *FT  A 
A9FTC  ■  SMALL A ( 9  »L 1 +SMALLA (  11 *L)*XI  *2.*SMALLA( 1?.L)*ETA 

A3ETC  «  SMALL A ( 3*L)+SMALLA(8»L)+XI  *(SMALLA(6.L 1 ♦ SMALL A ( 10.L 1*2. 1 
1  *ETA*  (  2.*  SMALLA  ( 6  ,L  1  -.SMALL A 1 11  .L  )  I 

SX(l.I)  *  SX(l.I)  +  C11*A2ETC  +  C12*A9ETC+C13*A3tTC 
SX( 3 . I )  *  SXI3.il  ♦  C33*A3ETC  +C13*A2FTC+C23«A9E TC 
120  SXI2.I)  *  SX(?.I)  +  1  ?*A2E T C  +  C22*A9ETC  ♦C23*A3ETC 

DO  130  J  *  1,3 

130  SX(J.I)  *  SX(J.I)  /  NOT 9 

r 

C  FIND  PRINCIPAL  STRESSES  PSTR 
C 


RHS  *  SORT!  (  SX(l.I)-  SX(2. !))••?  /4.  ♦  SXI3»I)**2) 

PSTR (1.1)  ■  I  SX(l.I)  ♦  SX(?,!M  /  2.  +  RHS 

PSTR  ( ? » I  )  «  (  SX(l.I)  +  SX(?,m  /  2.  -  RHS 

PSTR (  3  »  I  1  »9 .?  •*  ATANt?.*  SXI3.il/ISX  (l.I)-SX  (2.I111/PI 
200  CONTINUE 
C 

C  DERIVE  P  STAR  2 
C 

C  KS22INKD  X  N<D)  IS  IN  Sll.1,4) 

C  <12  (NUD  X  NKD 1  IS  ON  TAPF  10 

C  <21  IS  THE  TRANSPOSE  OF  <12 
C 
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N K?2  =  NX  D*NK  0 

CALL  ^XMUL  T  (SU  1 1 .4)  tFD<NUD+l  )  »PS2  tNKDt NKD .  1 ) 

REWIND  10 

DO  290  I  *  1 «M 

IF  <  I  .GT.l  I  GO  TO  23*' 

l DCFO  *  1 

K210IM  a  KRFM#NKD 

I  I  *  KRE* 

GO  TO  24^ 

230  II  *  N 

K21DIM  *  N*NKD 

LOCFD  *  KRFM+( I-?)»N+1 

74'<  continue 

READ  (PI  (0(  JJJ)  *  JJJ*1  ♦<71DTM) 

CALL  Tmxml'LIQ.  DEL  (LOCFD)  *PS2  *  1 1  tNKDt  1  ) 

29n  CONTINUF 
C 

C  PRINT  p-STAR  2 
C 

NPSTR  •  NKD/7 

IF(NKD.GT.NPSTR*7)  NPSTR  *  NPSTR+1 
DO  ^*>0  1*1  .NPSTR 

II  x  1*1 i-l )  +  l 

III  *  II  ♦  NUD 

12  =  MlNfM  I  1 +6  tNKD ) 

117*12+  NUD 

PRINT  1 03^ f  (  K •  K *  I  I  1 . 1  I  2  ) 

PRINT  1012«(PS2(<) t<* I  1 t  I  2  ) 

*60  CONT IN JE 

r 

C  PRINT  OUT  STRr 5SFS 

r 

I  =  1 

870  CONTINUE 
print  ioo:: 

L  INF  =  1 

871  CONTINUE 

DO  88^  J*1  *3 

T  F  (  ,)-?)  87*  .87^  ,P7* 

8  7**  PRINT  10d*  It  X(J).  SX(Jtl)  •  PSTRU.J) 

CO  TO  876 

874  PRINT  100?  •  Y  I  I  )  t  SX  (  J  .  I  )  ♦  PSTRUtl) 

GO  TO  876 

87*  PRINT  lOCRt  SX(Jtl)  .  PSTR(Jtl) 

876  CONTINUE 
8  8"  CONTINUF 

LINE  *  LINE  *  6 
I  =  1  +  1 
PRINT  10-74 

IF  (I  .GT.NN )  GO  TO  890 
IE  (LlNE.GT.bC )  GO  TO  870 
GO  TO  871 
po~  CONTINUF 
RE  T’  JRN 

1  on*  FORMATdHltflH  NODE  1 6X  t  8HPOS  I  T  1  ON  t  ox  ,  1 3HVECTOR  SIG^A.SXt 

1  1 8HPR I NC I °A L  STRESSES/) 

1  r'\  FORVATUH  7X«I3t6XtE17#^t2(6XEl*#8)  ) 

1  2  FORMAT  1 iH  3Xt3Xt6XtL12.b.2(6XElb.8) ) 

I  3  FORMAT ( 1h  3Xt3Xt6Xtl2X  .2C6XE15.8H 
1  ^  4  F  OR  M A  It/) 

ID] n  FOPMATt 1 H 1 .SoXt  13HDISPLACEMENTS#/ > 

1011  FORMATdH  *7  (  8X  •  4HDE  L  <  ?  I  3  •  1  H  )  )) 

I'M?  PORMATdH  t7(2X.El4.7>/) 

1030  FORMAT ( IH  71  8X «4HPS? < • I  3 1 IH)  )  ) 

END 
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SUBROUTINE  SMALLK(I) 


ANISOTROPIC 


this  subroutine  derives  the  matrix  small  k  for  the  triangle  i 

ALSO  CALCULATED  AND  STORED  FOR  THE  STRESS  ROUTINE  ARF  THE  A-STARS 

COMMON  tit  NO  (80*6).  FDI500),  NODE ( 500  . ? > «  X(?50).  Y(  ?50) . TH t ?50 ) , 

1  N.  KREM.  M.  S( 60*60*5 ) »  IDIR(?50.'t.  Ac.F(250).  T(12,12l, 

2  AMXU2.12).  BMX(  12*121 »  SP(6n.20).  ER(80),  RNIMRO). 

3  FNU(BO),  N0T<250.9).  C I J  ( 6 . 80  )  ,K  (  1 2 . 12 .80  )  ,  A$<6.6.80), 

A  G(6l>).  NKD*  NUD.  NN.  NT  . VF ( 80 ) . BE TA < 80  I 

REAL  KMX 

DIMENSION  C(3.3)  .  CPI3.3)  .  CDI3.3)  .  T1I3.H 
DATA  XLIMIT  /  •8E-B/ 

IA  ■  NO (1*1) 

IB  *  NO ( I  * 3 ) 

IC  ■  N0( 1.5) 

I AG  *  NOt I .4  ) 

IBG  *  NO ( 1.6) 

ICG  •  N0U.2) 

XIA  «  X ( I A  > 

X IB  «  X ( IB )  -  XIA 
XIC  «  XI IC)  -  XIA 
YIA  ■  Y ( I A  ) 

ETAB  ■  Y( IB )  -  YIA 
ETAC  »  Y(IC)  -  YIA 
XIO  »  (XIB+XICJ/3. 

ETAC  «  IETAB+ETAC)  /3. 

DELT  «  XIB*ETAC  -  XIC*ETAB 

T1 1 ■  -(TH( IA)»(ETAC-ETAB)  -  TH( IB)*ETAC  ♦  TH ( IC  l*FT AB  l  /  DFLT 
T2  ■  -  ( TH( I A ) * ( X  IB  -  XIC)  *  TH( IB ) *X IC  -  TH(TCI*XIB  )  /  OFLT 
TO  «  TH(IA) 

EM2  ■  FTAC  /XIC 
IFIXIB.NE.O. )  GO  TO  15 
EMI  «  EM2 
GO  TO  17 

15  EMI  «  ETAB  /  XIB 

17  EM3  «  (FTAC  -  ETAB)  /  (XIC  -  XIA) 

B3  -  ETAB  -  EM3  *  XIB 

B32  *  B3**2 

B33  *  B3?*B3 

B3A  *  B33*B3 

EM? 1 1  •  FM2  -  EMI 

FM212  *  (EM2+EM1)  *  EM211 

EM? 13  *  EM?**3  -  EM1**3 

EM23)  «  EM?  -  FM3 

EM?1A  «  FM2#*A  -  EM1**A 

X IC3  *  XIC  «*  3 

X  ICA  «  XIC  *  X IC3 

XIC5  ■  XIC  *  XICA 

XIB3  *  XIB  **  3 

XIBA  «  XIB  *  X IB3 

X IB5  »  XIB  •  XIBA 

FMFC1*  FM3*XIC  +  R3 

FMFB1 *  EM3»XIB  +R3 

EMECB3  «  FMEC 1**3 

EMECBA  «  EMEC1  *  EMECB3 

EMECB5  *  EMEC1  •  EMtCBA 
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EMER33  *  EMEB 1**3 
FMEBB4  «  FMEB1  *  EMEBB3 
FMEBB5  ■  EMEB1  *  EMEBB4 
FCR  1  *  XTC-XIR 
FCB?  •  X l C**2  -  XfR**? 

ECB  3  *  XIC**3  -  XIB»»3 
ECB4  =  XIC**4  -  X I  B#*4 
FCB5  *  X  IC## 5  -  XIB**5 
EM32  *  EM3  *  EM3 

EM3 3  ■  EM3  *  EM32 

EM22  •  EM2  *  EM2 

FM23  •  EM2  *  EM2? 

EM24  *  FM?  *  EM23 

C 

C  CALCULATF  THF  GAMMAS 

C 

GO  «  DFi_T  /  2* 

G1  *  GO  *  XU: 

G2  ■  GO  *  ET AO 

G3  «  FM?U  #  XIB4/4.  ♦  FM231/4.#FCB4-B3/3.  *  ECB3 

G6  ■  EMp  1 1  *  XIBS/S.  ♦  FM231/«>.*ECB5  -  B3/4.  *  ECB4 

C 

IF C FM3.NE.G. )  GO  TO  18 

G4  *  EM213  *  X  I B4/1 2 • ♦  EM23/ 12 .*ECB4-B3 3*ECBl/3 . 

C 

G5  *  EM212  «  XIB4/8.  «■  EM22/8.  *ECB4-B32*ECB2/4 . 

C 

G7  •  FM212  *  XIBV10.+  FM?2 / 1 0.#FC*3-B32*FCB3/6 • 

C 

G8  ■  EM213  *  XIB5/15.+  FM23 / 1 3 •  *ECB5-B3?*ECB2/6 • 

C 

G9  * ( EM2 14*X I 85  ♦EM24*ECB5 ) /20.  -B34*ECBl/4. 

C 

GO  TO  19 

r 

18  G4  *  FM?  1 3  #  XIB4/12.+  EM?V12.*FCB4-(FmECB4-FMFBB4>/(12.*FM3) 

GS  «  FM?  1  ?  *  XIB4/P.  ♦  FM22/8,  #FCB4~  ( rMErB4-FMFBB4 )  /  ( 8 .*tM3?  > 

1  ♦  03/ ( 6.*FM32 )  *  ( EMFCB3  -  EMEBB3 ) 

C 

C 

G7  -  EM212  *  X  I B5/10.+  EM22/10.*ECB5- ( ( EMECB5-EMEBB5 ) /10. 

1  -  ( EMECB4-EMEBB4 ) *H3/4 •  ♦  B32* ( EMECB3-EME8B3 ) /6. ) /EM33 

r 

Gfl  =  Ey?  1  3  *  XIBS/15.+  EM23/1S#*FCBS-  ( EMFCBS-EMFPB5 )  / ( 1 S  **EM32 ) 

1  +  B 3 / ( 12* #EM3? )  *  IEMFCB4-FMEBP4) 

C 

G9  «(FM214#XIB5  +EM24*FCB5  -  ( FMFCBS-E ME BBS )  /EM3 )  /  20. 

C 

r 

19  nun  =  cos  calf  <  m  ) 

T  ( ?  ♦  ? )  *  COB ( Al F (  JR) ) 

T  (  3  *3  )  *  COSCALFC I C )  ) 

T  (  4  1 4 )  *  COSCALFC  I  AG  I  ) 

T  C  S  ♦  5 )  *  COSCALFC  1BG)  ) 

T(6.6)  *  COSCALFC  ICG )  ) 

TC7.7)  =  T( 1 *1  ) 

T (8*8)  *  T  (  2  ♦  2  ) 

T<9*9)  =  T  <  3  t  3  ) 

T(1!M0)  *  T  (4.4) 

T(lUll)  «  TCS.5) 

T C 1 2 *  1 ? I  *  TC6.6) 

T  ( 1  ♦  7 )  *  -SINC ALFC IA  )  ) 
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T(2,8)  -  -SINIALFUB  )  ) 

T (3 .9 )  •  -SINIALFIIC  )  ) 

1(4.10)  -  -SINIALF(IAG)  ) 

1(5.11  )  «  -SIN(A(_F(  IBG)  ) 

TI6.12)  »  -SINIALFI ICG)  ) 

T ( 7 . 1 )  *  -  T( 1.7) 

T  ( 8  .  ? )  ■  -  T  (  2 . 8  ) 

1(9.3)  *  -  TI3.9) 

7(10.4)  -  -  TI4.10) 

Till. 5)  *  -  TI5.11) 

7(12.6)  «  -  TI6.12) 

VR  «  1.  -  VF( I) 

GR  «  FR ( I ) /(2.*< l.+RNUI I ) )  ) 

GF  ■  EF ( I ) /(2.*( l.+FNUI  I)  )  ) 

EXP  ■  ER  < I )*VR+EF( I )*VF (  I ) 

EYP  -  ER(  I  )*EF(  I  )  /  ( b  R ( I )*VF  ( I ) +EF ( I ) *VR ) 
XYNUP  ■  RNIJI  I  )  *VR  +  FNU (  I  ) *VF  I  I  ) 

GXYP  »  GR*GF  /  ( GR*VF  (  I  )  +  GF*VR  ) 

B11P  «  1 • /EXP 

B12P  «  -XYNUP*B11P 

B2?P  «  l./FYP 

B33P  «  l./GXYP 

DELTA  ■  B11P*B22P-B12P*B12P 

CPI  1.1)  «  B22P/  DELTA 

CPI  2  » 1 )  *  -  B12P/  DELTA 

CPI3.1)  ■  0. 

CPU. 2)  *  CPI2.1) 

CPI  2.2 )  -  B11P/DFLTA 

CPI3.2)  ■  0. 

CPU. 3)  ■  0. 

CPI ? . 3 )  «  0. 

CPI  3.3 )  -  1./B33P 
EL  «  COS(BETA( I ) ) 

EM  ■  SINIBETAI  I  ) ) 

Tl(l.l)  ■  EL**2 
TK1.2)  ■  EM**  2 
T 1  ( 1.3)  *  FL  *  EM 
TK2.1)  ■  Till. 2) 

T1I2.2)  *  Tl(l.l) 

T1I2.3)  »  -Till. 3) 

T1I3.1)  ■  2.*T1 (2.3) 

Til  3.2)  «  -TH3.il 
TK3.3)  »  TK1.1)  -  Till. 2) 

T2  NOW  IN  T 1 

CALL  MXMULT (CP.T1.CD.3.3.3) 

T1I3.2)  -  T1 12 .3 ) 

TK3.1)  »  Till. 3) 

Till. 3)  ■  2.*T1 ( 3.2 ) 

T 1 ( 2 * 3 )  »  2.*T 1(3,1) 

T1  INVERSE  NOW  IN  T1 

CALL  MXMULT (T1.CD«C,3,3«3) 

NOW  HAVE  C  MATRIX  IN  C 


CIJ(l.I)  -  C(l.l) 
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n  o  o  o  ^ 


c i j ( ? «  m  *  c  ( 1 .2 ) 
CIJMiI  )  =  C  (  1*3) 
C  I J  (  4  *  I  )  *  C  (  ?  «2  1 
CUM. I  )  =  C<2,3) 
UJ(6tl)  *  C<  3*3) 


CALCULATE  A-STAR  MATRIX 


A  S  {  1  ♦  ?  •  I  )  * 
AS(  1  •  3  *  I  )  = 
AS( ?♦?,!>  « 
AS ( 3 , 2 ,  I  )  * 
ASIA. 2*1 )  * 
AS ( 5*2* I )  c 
AS ( 6 *  2  *  I  )  * 
A  S  (  2 . 3  *  I  )  = 
A  S  (  3  *  *  ,  T  >  = 

AS  < 4.3. I  )  * 
ASM.  3,  I)  = 
AS ( 6 , 3 , 1  )  ' 
DO  20  J  = 
AS ( J.l . 1 )  * 

AS ( J  *  4 . I  )  * 

A  S  (  J  .  5  .  I)  = 
2^  AS( J.S.I )  * 


0  • 

D. 

X  IR 

X  IC 

XI  I  AG )  -  XIA 
X < IBG)  -  XIA 
X ( ICO)  -  XIA 
ETAR 

F  T  AC 

Y (  I  AG)  -  VIA 
Y ( I BO )  -  t  i  « 

Y(  ICG)  -  YIA 
♦  6 
1. 

A  S ( J  *2 .  I  )**2 
AS  < J .2 . I )  * 
AS(J»^  I  )**? 


AS( J.3.  I) 


CALL  INVERT!  AS(  1.1.  I  ) .6,12*XLIMIT.FLAG) 
IEtELAG.NF.r • ) GO  TO  SOO 


NOW  HAVE  A-  STAR-INVERSE  .  GET  A-MATRIX 
DO  25  J*  1  ,S 

on  2*"  K=  1  .f 

AMX(J+6.K*6)  =  AS(J.K.I) 

25  AMX ( J  ,<  )  *  A S ( J  * K  .  I  ) 

TGI  * T 1 1  *  G1+T2*02+T(*G0 

TG3  ’Til*  G3+T?*G5+T  >*G1 

T G5  s  T  1  1  *  G5  +  T2*G4  +  U'*G2 

TG6  =  T  1  1  *  G6+T?*G7+T >*G3 

TG7  ’Til*  07  +  T  ?*Gfl*TO*G5 

TOR  =  T  1  1  ♦  G  R  +  T  2 * G 0 -f  T 0 * G4 

R MX (  2.7 )  *  C (  1  .1  )  *  TGI 

RMX (3.2)  *  C<1  *  3 ) *  T  G 1 

HMX (4.2 )  ’  2  •  *  C ( 1 . 1  )  *  TG3 

0Mx  (  5  ,  ?  )  *  CU.l)  *  TGr>  ♦  C  (  1  *  3  )  *TG3 

RMX  (6.2  )  =  2.*C(  1 .3  )  *TG*' 

RMX ( R , 7 )  =  PM.xn,?) 

RVX (4,2)  *  C ( 1 . 2 )  *  TGI 

°MX  (10.2)  *  Pm  *  CM, 3)  ♦  TG* 

HMX (  1 1  , 2 ) =  C (  1 , 2 )  *  T  03  +C (1,3)  *  TG5 
RMx ( 12.2 ) *  ?•*  C ( 1 »?  )  *  TG5 
P WX ( 3 , 3 )  *  C (  3,3)  *  TGI 

HMX (4,3 )  =  RMX ( 10,2 ) 

B  M  X  (  5 , 3  )  =  C (3,3)  *  TG3  +C  ( 1 , 3 )  *  TG5 
tjvx(6M)  *  2**  C  (  3 , 3 ) *TG5 
HMX(R.I)  =  HMX (3,3 ) 

°vx (9,3)  =  f ( 2 ,3 )*TG) 

RMx (  1,3)  =  2 • *C (3*3) #T  G3 

HMxnit’l*  C  (  3 , 3  )  ♦  TG5  ♦  C(?,3)  *  TG  3 
RMXM2.3)  =  2.*C(2,3  )*Tr-S 
9*X (  4,  4)  =  4**  C (  1,1)  *  TG6 

RMX  (  5,  4)  =  ?.*  C(l,n  *  TG7  +2.*C  (1  ,3)*TG6 
RMX (5,4)  =  4 • *  C (  1,3) *  TG7 

°VX(R,4)  =  RMX ( 1 C ,2 ) 

Rvx (  9,  4)  *  2.*  C (  1  , 7  )  *  TG3 

H VX (  ] 0 , 4  )  *  4**C(1 ,3 )*TG6 

RMX(11,  4)  x  C(l. 2)  *  TG6  ♦  2 • *C ( 1 ♦ 3  )  *TC7 
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RMXI12*  4)  *  4.«  CI1.2)  *  TG7 

BMX  (  5.  '»'»  *  C  ( 1  •  1  J*TG8  +  CO. 3)  *  TG6  ♦  2.  «■  CI1.3)  *  TG7 
BMX(  6.  «.  i  *  2.*  C(  3 .3  )*TG7  2.*Cfl.3)  *  TC.B 

BMX (  8.  5)  «  BMX (5.3) 

BMX (  9.  5)  •  C  ( 1 .2 )*TG5  +  CI2.3)  *  TG3 

RMXIin.  5 )  *  ?.*  CO. 3)  *  TG6  +  2.*C( 1 .3)*T07 

BMX  (11.  6 )  =  (C(l>2)+C(3.3)  )*TG7  +C ( 1 . 3 ) «TG8*C ( 2 . 3 )#TG6 

9  MX  (  12«  5)  *  ,?•*  C(1  *?  I*  TGB  +  2.*C  (  2. 3  )  *TG7 

BMX  (6.61  =  4.*C( ?.3>*TG8 

BMX  (  8.  6)  «  Bj<lX (-6.3) 

BMXI9.6)  *  BMX  f  1 2  *  3  ) 

9MXI1J.6)  *  4  »*C ( 3 .3 ) *TG7 

BMX (11.  6)  «  2.*  C  (  3 . 3 ) *  TG8  +  2.*C I  2 .3 ) *TG7 


BMX (12.6)  ■  4.*C(2.3 )#T08 
BMX (  8.  B)  *  BMXIB.3) 

BMX ( 9 . 8  )  =  BMX (9.3) 

BMX  (  1  .^ »  8)  =  BMXI10.3) 

BMX (11.  8)  =  BMX (11.3) 

BMX  (12.8)  =  BMX  (12  0) 

BMX (  9.  9)  »  C ( 2 .2 ) *TG1 


BMX(1'),9  )  =  2.*C ( 2  *  3 ) *TG3 
BMX (11.  «)  *  C ( 2 .2 ) *TG3  +  C(2.3)*TG5 
BMX(ll.lO)  *  2  .*C  (  3.3  )*TG7  ♦  2 . *C ( 2 *3 >* TG6 
BMX ( 12.9)  *  C ! 2  »  2 ) *TG5  *2. 

BMX  (  11.  lx)  *  COO)  *  TG6  +  C(3.3)  *  TG8+2  .*C  ( 2 .3 )  *TG7 
BMX ( 1 0  .  ID )  *  4.*C ( 3  »3 )*TG6 
BMX (12.10)  «  4»*C (2.3) *TG7 

BMX  (  12*  11)  ■  2.*  C ( 2  *2  )  *  TG7  2 .*C (  2 . 3  ) *TG8 

BMX (  12.  12)  >  4.*  C (2.2 )  *  TGB 


APPLY  SYMMFTRY  HERE 

DO  200  J*3  *  1 2 
J1  ■  J-l 
DO  200  L*2* J1 
200  BMX(L.J)  «  BMX(J.L) 


CALL  MXM!)LT(B«X.AMX.XMX(1.1.I  +  1)  .12.12  ,  12) 

CALL  MX MOLT (KMXI1.1.I+1  )  »T  .KMX (1*1* I). 12. 12. 12) 

DO  25r  J=1  .6 
DO  250  K-1.6 

250  AMX(J.K)  =>  AMX  (  K+6  * J+6  ) 

DO  260  J=1 .6 
DC  26''  K*1.6 
AS(K.J.I)  «  AMX(J.X) 

260  AMXIJ+6.K+6)  *  AMX(J.K) 

TRANSPOSE  OF  INVERSE  OF  A  NOW  IN  AMX.  INVERSE  OF  A-STAR  NOW  IN  AS 


DO  270  J  =  1  .6 
T ( J  *  J+6  )  *  — T(J.J+6) 

270  T ( J+6 ♦ J  )  =  -TIJ+6.J) 

INVFRSF  OF  T  NOW  IN  T 

CALL  MXMULT ( AMX  .KMX (1.1.1  ) .KMX ( 1  .  1 . I  + 1  )  ♦  12 . 1 2 . 1 2 ) 

CALL  MXMIJLT(  T.KMxtl.l.I  +  1  )  .  KMX  (  1  .  1  .  I  )  .  1  2 . 1 2 . 12  ) 

C 

C  K-PRIME  MATRIX  NOW  IN  KMX(l.l-I) 

C 

RETURN 

500  PRINT  100'. I 

1000  FORMAT ( 44H1  COULD  MOT  INVERT  A-STAR  MATRIX.  TRIANGLF  .13) 
STOP 
END 
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SUBROUTINE  INVtRT(B.K.K2.XMIN.ELAG) 


r  this  suproutinf  sfts  up  a  unit  matrix  adjacent  to  bik.k) 

C  ELEMENTARY  ROW  OPERATIONS  ARF  THEN  PERFORMED  ON  THE  NEW  X  X  2K  MATRIX 
C  Tn  REOUCF  fl(K.K)  TO  A  UNIT  MATRIX.  THIS  WILL  PLACE  THE  INVERSE  OF 
C  The  MATRIX  BU.K)  IN  Tn>  RIGHT  half  OF  B(K.2K> 

C  ON  EXIT.  THE  INVERSE  OF  H  REPLACES  B 

C  19+S  N  ARRAY  OF  2*K**2  LOCATIONS  CONTAINING  THE  MATRIX 
f  X  ISTHF  D I  M*"NS  I  ON  OF  THE  SQUARE  MATRIX  A 
C  K?  IS  ?«K 

r  XV I N  IS  THE  SMALLFST  ALLOWABLE  MAGNITUDF  OF  THE  PIVOT 
C  FLAG  WILL  BE  RETURNED  AS  0.  IF  THE  INVERSION  WENT  OFF  OX 
C  FLAG  WILL  PF  RETURNED  AS  10.  IF  A  PIVOT  ELEMENT  WAS  TOO  SMALL 
C  FLAG  SHOULD  BE  TESTFO  AFTER  EACH  CALL  TO  THIS  ROUTINE 

C 

DIMENSION  BIK.X2) 

C 

FLAG  *  0. 

C 

C  SET  UP  UNIT  MATRIX 

f 

IF(X.GT.l)  GO  TO  20 

IFIABSIBI 1.1 ) ) .LT.XMIN)  GO  JO  10 

B  C 1  •  1  >  »  l./B< 1.1) 

RETURN 
2"  CONTINUE 
DO  1  J ■ 1 «K 
DO  1  J»1.X 

B ( I  .X>J)  ■  ft. 

I F ( I.EO.J)  B I  I .K*J)  ■  1. 

1  CONTINUE 
C 

C  FIND  LEADING  ELEMENT  WITH  GRFATFST  MAGNITUDF 

r 

DO  A  J-l.X 

M  u  J 
N  ■  J+l 

IFU.EO.XI  GO  TO  21 
DO  2  L*N.K 

IF  (ABSIB(M.J) ) .LT.AB SIBIL. J) ) )  M«L 
?  CONTINUE 
?1  CONTINUE 

IF  (ARSCniM.J) ). LT.XMIN)  GO  TO  10 
IFtJ.EQ.X)  GO  TO  31 

C 

C  INTERCHANGE  JTH  AND  MTH  ROWS 
C 

DO  3  L-J.K2 
D  »  RIJ.L) 

B(J.L)  *  B(m,L) 

B(v,l)  »  D 

3  CONTINUE 
31  CONTINUE 

C 

C  2FRO  OUT  PIVOTAL  JTH  COLUMN.  SKIPPING  PIVOTAL  JTH  ELEMFNT 

C 

C  DtVIDE  JTH  POW  BY  PIVOT 

r 

no  4  m.n.x? 

Bi J.M)  «  H(J.M)  /  B ( J . J ) 

4  CONTINUE 

DO  6  m*i ,< 
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n  n  o  ri  o 


c  M  DETERMINES  ROW  BEING  MOD  IF  I  ED t  ONE  WHOLE  ROW  AT  A  TIME 
C 

IF  (  M.EO.J  )  GO  TO  f, 

DO  •>  L  *'•»<? 

C 

C  L  DETERMINES  ELEMENT  IN  THE  MTH  ROW 
C 

B(M*L)  ■  B(M*L)  -  BIM.J)  *  R(J*L5 

5  CONTINUE 

6  CONTINUE 

r 

C  INVERSE  OF  B  IS  NOW  IN  RIGHT  HALF  OF  R(K,K2) 

C  NOW  MOVE  B  INVERSE  TO  WHERE  B  WAS 
DO  7  I *1 »K 
DO  7  J«1»K 

B ( I ♦ J )  ■  B ( I »J+K) 

7  CONTINUE 
RETURN 

in  FLAG  =  10. 

RETURN 

END 


SUBROUTINE  MXCON < A,B »X *M ,N ) 

C 

c  this  surroutinf  multiplies  matrix  a  imxni  by  constant  x,  resjlt  in  r 
C  A  MAY  be  sa^f  as  b. 

C  THIS  SURROUTINF  multiplies  MATRIX  a  (MXN)  by  constant  x,  result  in  b 
C  A  MAY  HE  same  AS  3. 

DIMENSION  A(M.N)  *  B ( M»N ) 

DO  1  I-l.M 
DO  1  J-l.N 
B ( I  * J ) *  X *A ( I .  J) 

1  CONTINUE 

return 

FND 


SUBROUTINE  MXSUB(A,B.C»M.N) 

this  subroutine  subtracts  matrix  b  from  matrix  a. stores  result  in  c 

A,  B.  AND  C  ARE  (M  X  Nl  (C  CAN  BE  THE  SAME  AS  A  OR  B) 

DIMENSION  A(M*N)  »  B(M»N)  tC(M»N) 

C 

DO  1  I  *  1 »M 
.DO  1  J* 1 »N 

C(I.J)  «  A ( I  * J )  -  BIT  « J) 

1  CONTINUE 
RETURN 
END 
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4 


SURROUT  INF  m<vui.T  (  A  «B  «C  «N’*N  «X  ) 

(  this  SlIRPOUTlNF  VUL  T I  PI  ! r  S  MATRIX  A  BY  matrix  «  *NI'  STOOFS  THF 

r  PRODUCT  IN  C.  (C  CANNOT  BE  THF  SAME  AS  A  OR  B.) 

C 

C  A  IS  (M  x  N) 

C  R  IS  ( N  X  X ) 

C  C  IS  (mx  X ) 

c 

DT^FNMPN  A(v,N)  •  Q(\*K)  « 

r 

no  l  i«i»m 

00  1  L*lt< 

C  l  I  *L  )  =  • 

DO  1  J  =  1  *  N 

C  (  !  *  L  )  *  C(hL)  ♦  A  (  !  t  J  )  * 

1  CONTI  NOT 
RETURN 
F  NO 


SURROI  r I Mt  TMXMUL (A'R'C'NtM'K ) 

r 

f  THIS  $W*ROOTINc  Mt*L  T I PL  I r  s  matrix  *  *Y  THP  transpose  of  matrix  a 
r  thf  PRoourT  is  ahofp  r 

r 

r  a  is  in  x  mi 

C  H  IS  (N  X  X) 

C  CIS  I m  x  x I 

c 

Dive  NS  ION  A(N.M)  ,  B(N,X)  «  C(V.X) 

Cp  1  !=1. v 

'AO  1  l  =1  ,X 
r  C  (  T  *  L  )  »  . 

DO  1  J- '  .N 

C(I.L)  =  C I  I  *  L  )  +  A(J.I)  *  P(J.L) 

1  CONTINlJF 
RE  TURN 
FND 
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